In this study, we investigate the implementation of a Proper Orthogonal Decomposition (POD) Polynomial Chaos Expansion (PCE) POD-PCE surrogate model for the propagation and quantification of the uncertainty in hydraulic modelling. The considered model consists of a system of multilayer shallow water equations with a mass exchange between the layers and over stochastic beds. As a numerical solver, we propose a finite volume characteristics method that does not require eigenstructure of the system in its implementation. The method is fast, accurate and can be used for both slowly and rapidly hydraulic simulations. The propagation and influence of several uncertainty parameters are quantified in the considered numerical methods for multilayer shallow water flows. To reduce the required number of samples for uncertainty quantification, we combine the proper orthogonal decomposition method with the polynomial Chaos expansions for efficient uncertainty quantification of complex hydraulic problems with a large number of random variables. Numerical results are shown for several test examples including a dam-break problem over a flat bed, and a wind-driven recirculation flow on flat and non-flat bottoms. Results are also presented for the case study of a recirculation flow problem in the Strait of Gibraltar. The results demonstrate the robustness of the uncertainty quantification method compared to the standard Monte-Carlo simulations. The results presented in this study suggest that the use of surrogate modelling may save a considerable amount of the necessary computational cost for all the considered cases.