Tackling random fields non-linearities with unsupervised clustering of polynomial chaos expansion in latent space: application to global sensitivity analysis of river flooding

A surrogate model is developed to accurately approximate a two-dimensional hydrodynamics numerical solver in order to conduct a reduced-cost variance-based global sensitivity analysis of the hydraulic state. The impact of uncertainties in river bottom friction and boundary conditions on the simulated water depth is analyzed for quasi-unsteady flows. An autoencoder technique adapted to non-linear variable dimension reduction is used to reduce the multi-dimensional model output so that the formulation of the surrogate remains computationally parsimonious. In addition, following the divide-and-conquer principle, a mixture of local polynomial chaos expansions is proposed to deal with non-linearity in the hydraulic state with respect to uncertain inputs. Machine learning techniques are used to automatically partition the input space into clusters that are not affected by non-linearities and support accurate surrogates. This combined strategy is applied to a reach of the Garonne River where river and floodplains dynamics are simulated by the numerical solver Telemac-2D. The merits of this strategy are highlighted when the flood front reaches regions where the topography features a strong gradient and where, consequently, strong non-linearities occur between the water depth and friction as well as hydrologic input forcing. By applying this strategy, the Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q_2$$\end{document} metric improves by 90% compared to a classical polynomial chaos expansion surrogate, resulting in a much more reliable sensitivity analysis. This is particularly important in floodplain areas where human and economic activities are at stake.


Flood monitoring
According to the World Health Organization (WHO), in Europe, floods are the most common natural hazard leading to emergencies, causing extensive damage, disruption and health effects (WHO 2017). Over the last 20 years, flood events have been recorded in 49 of the 53 member states. Estimates from WHO Regional Office for Europe, based on data from the international disaster database (EM-DAT), indicate that approximately 400 floods have caused the deaths of more than 2000 people, affected 8.7 millions others, and generated a loss of at least 72 billion euros over 2000-2014(Guha-Sapir et al. 2015. The magnitude of the physical and human costs of such events can be reduced if adequate emergency prevention, preparedness, response, and recovery measures are implemented in a sustainable and timely manner (WMO 2013). Resilient and proactive health systems that anticipate needs and challenges are more likely to reduce risks and respond effectively during emergencies, thereby saving lives and alleviating human suffering. In this sense, several measures have been taken by governments and environmental organizations to minimize these effects (EFAS 2017), including the assessment and mapping of flood and tsunami health risks in order to indicate areas at highest risk, identify and analyze capacities for flood risk prevention, preparedness, response, and recovery with respect to the assessed flood risk, determine recommended actions for flood health emergency risk management, and assess resources and identify priorities for action.
The climate community estimates that about 1.3 billion people will be affected by flooding by 2050 due to climate change, increase in population density, and global degradation of environmental conditions (Arnell and Gosling 2016). It is increasingly clear that climate change has detectably influenced several water-related variables that contribute to floods, such as rainfall and snow melt. As global warming contributes to exacerbating sea level rise and extreme weather, floods are expected to grow by approximately 45% by the end of this century (Kulp and Strauss 2019). Thus, it is crucial to understand, assess, and anticipate flood events.
Flood monitoring benefits from world wide efforts by international programs dedicated to Earth observation from space, such as Copernicus, as well as from space agencies that support missions, such as Sentinel, or Surface Water Ocean Topography (SWOT) designed to study the topography of oceans and continental bodies of water (Biancamaria et al. 2016). In spite of the increasing volume, resolution, and precision of remote sensing water surface elevation observations, the prediction of flood events requires the use of reliable and robust numerical hydrodynamic models.
In France, the forecasting and vigilance of hydrological events likely to generate floods is ensured by Service de Prévision des Crues (SPC) whose action is coordinated by Service central d'hydrométéorologie et d'appui à la prévision des inondations (SCHAPI) of the Ministry of the Ecological Transition. The SPC/SCHAPI network works in partnership with Météo-France, which provides it with the meteorological variables (observations and forecasts) necessary to drive their hydrodynamic models.

Hydrodynamic numerical solvers
River hydrodynamic models are used to predict river water depth and discharge from which flood risk can be assessed. These predictions provide a Decision Support System (DSS) (Daupras et al. 2015) with informed hydraulic parameters and variables (water depth, discharge, and velocity) along with their evolution in the future for leadtimes that range from a couple of hours to a couple of days depending on the dynamics of the catchment. DSS are thus able to manage flood risk and eventually issue alerts for protective actions. Several research projects and concerted actions have been funded on the subject of river flood monitoring. For instance, the Hydrologic Ensemble Prediction EXperiment (HEPEX) aims to develop and demonstrate new hydrologic forecasting technologies and to facilitate the implementation of beneficial technologies into the operational environment (Schaake et al. 2006). The European Flood Awareness System (EFAS), initiated in 2003 (Thielen et al. 2009), seeks to improve flood preparedness in transnational European river basins by providing medium-range deterministic and probabilistic flood forecasting information, from 3 to 10 days in advance, to national hydro-meteorological services, e.g., SPC and SCHAPI.
Hydrodynamic numerical models are generally based on a deterministic approach that solves the Shallow Water Equations (SWE) derived from the free surface Navier-Stokes equations (de Saint-Venant 1871;Sohr 2001) and are prone to uncertainties. The uncertainty in the water depth and discharge field computed with a hydrodynamic solver is due to uncertainty in simplifying assumptions with respect to physics, particularly with respect to the flow dimension, approximate knowledge of hydraulic parameters, and imperfect description of forcing and geographical data. Uncertainty quantification aims to quantify and rank the major sources of uncertainties, thus allowing for a better informed and, eventually, improved hydraulic forecast.

Surrogate models for sensitivity analysis
Global Sensitivity Analysis (GSA) consists in studying how the uncertainty in the output of a model (numerical or otherwise) can be apportioned to the different sources of uncertainty in the model input (Saltelli 2002;Razavi et al. 2021). The aim of GSA is to identify and rank the parameters that contribute mostly to the variability of the output of a model, also called a Quantity of Interest (QoI). It thus identifies which source of uncertainty should be reduced to most efficiently reduce uncertainty in the simulated QoI. A popular approach for sensitivity analysis is based on the decomposition of the output variance as the sum of the contributions associated with each input parameter and their combinations from which Sobol sensitivity indices are computed (Archer et al. 1997;Saltelli 2010). Extensions of those indices exist in the case of functional output (De Lozzo and Marrel 2017). This approach thus relies on sampling the uncertainties in the input space and the propagation of uncertainties through the model. Monte Carlo (MC) simulation is the most common technique used for sampling and Sobol indices computation (Sobol' 2001). However, its convergence is slow as it scales inversely to the square root of the MC sample size and its cost becomes prohibitive for computationally expensive models such as two-dimensional (2D) hydrodynamic solvers, especially in the context of realtime forecasting. To overcome this limitation, surrogate models may be used in place of the direct solver (Razavi et al. 2012). A surrogate model is a cheap-to-evaluate and parsimonious data-driven emulator of a reference model. This reference model can be seen as a black box that only provides a limited number of evaluations or observations. Thus, its output is known only at a few selected input points by means of a design of experiments. Then, the surrogate model seeks to approximate the reference model from this sparse input-output dataset. A variety of approximation techniques have been developed and applied as surrogates, such as linear regression models (Haldar and Mahadevan 1999), multidimensional scaling (Kruskal and Wish 1978), splines (Friedman 1991), Gaussian process (Rasmussen and Williams 2006), radial basis functions (Buhmann 2003), polynomial chaos expansions (Ghanem and Spanos 1991;El Garroussi et al. 2019), and artificial neural networks (Kasiviswanathan and Sudheer 2013). Some of them can interpolate the learning input-output dataset, e.g., Gaussian process regression, whereas others are designed to model the relationship between a QoI and sources of random uncertainty, e.g., polynomial chaos expansions.
The surrogate model based on Polynomial Chaos Expansion (PCE) (Lucor et al. 2004;Le Maître and Kino 2010) has proven useful in a wide range of applications, providing a low-cost yet accurate meta-model to estimate sensitivity indices (Sudret 2008;Crestaux et al. 2009). This surrogate model relies on the decomposition of the output random variable onto an orthonormal basis of polynomial functions. The polynomial coefficients are obtained either by using intrusive methods requiring access to the analytical code behind the numerical solver (e.g., Galerkin projection) or non-intrusive methods that rely on a learning database using the numerical solver as a black box (e.g., least square approximation). For steady flow in 1D and 2D, Roy et al. (2018), Goutal et al. (2018), andEl Garroussi et al. (2020) show that the PCE surrogate model succeeds in representing the response in water depth to uncertainties in river bottom friction and upstream discharge, allowing for an efficient computation of Sobol indices, water depth Probability Density Function (PDF), and water depth error covariance matrix over a reach of the Garonne River in southwest France.
However, PCE surrogates tend to struggle when applied to problems that feature non-polynomial non-linearities (Li and Ghanem 1998) or stochastic discontinuities that may occur for time-varying processes (Najm 2009). Indeed, for unsteady flow with a 2D hydrodynamic model, strong non-linearities in the water depth response to changes in bottom friction and upstream discharge may occur when water overflows the minor bed of the river; especially near dikes and in areas where bathymetry features strong spatial gradients. These non-linearities tend to exacerbate in unsteady regime, when the flood front, characterized by a non-zero velocity and a zero water depth, enters a previously dry floodplain domain. In this context, classical PCE meta-modeling is no longer adequate El Garroussi et al. 2020). Different approaches with varying degrees of complexity have been proposed in the literature to address the issue of PCE meta-modeling in the presence of non-linearities. Examples include multi-resolution/multi-element polynomial chaos expansions (Le Wan and Karniadakis 2005), regression trees (Torre et al. 2019;Choubin et al. 2019;Marelli et al. 2021), multivariate adaptive regression splines (Friedman 1991;Dertimanis et al. 2018), among others. They rely on the idea of partitioning the input parameter space into (often disjoint) sub-spaces followed by the use of intrusive or non-intrusive methods to estimate PCE coefficients. The surrogate model strategy should also be compatible with the dimension of the numerical solver output. Indeed, for functional output discretized over a mesh grid, the construction of a surrogate per mesh node would be computationally expensive, and could potentially lead to inconsistency as spatial coherence of the signal simulated field is not accounted for. The dimension of the model output should thus be reduced before the meta-modeling algorithm is applied (Bellman and Kalaba 1961;Lataniotis et al. 2020;El Garroussi et al. 2019). Dimension reduction stands in the transformation of high-dimension data into a meaningful representation of reduced dimension. On one hand, linear strategies, such as Principal Component Analysis (PCA) (Wold et al. 1987), linear discriminant analysis (Izenman 2008), factor analysis (Yong and Pearce 2013), and 3-way tables (Cichocki et al. 2009) are often used. On the other hand, kernel PCA (Schölkopf et al. 1997), Laplacian eigenmap (Belkin and Niyogi 2003), locally linear embedding (Roweis and Saul 2000), isomap (Tenenbaum et al. 2000), and AutoEncoder (AE) (Wang et al. 2016) are used to deal with non-linearities within data.

Objective and outline
In this paper, a surrogate model is developed to represent the 2D water depth field over the river and floodplain of the Garonne river, with respect to bottom friction and discharge. The surrogate model strategy aims to overcome the limitations of the classical PCE approach from El Garroussi et al. (2020), which provides a poorly predictive surrogate model in the presence of non-linearities for a transient flow. Both PCA and AE algorithms are investigated to reduce the dimension of the hydraulic output field so that the computational cost of the surrogate construction remains parsimonious. A Mixture of Polynomial Chaos Expansion (MPCE) approach is then implemented in the reduced space. Machine Learning (ML) techniques are used to partition the input space into disjoint clusters that are not affected by non-linearities and support an accurate PCE surrogate. The overall strategy, further denoted as reduced Mixture of Polynomial Chaos Expansions (rMPCE), allows to take advantage of the advances made in PCE surrogate modeling for local regression as well as in ML for dimension reduction and clustering. The resulting surrogate is used to carry out a GSA in order to rank the sources of uncertainty with a variance-based sensitivity analysis in the presence of non-linearities and at a parsimonious computational cost. The rMPCE approach and its application for the computation of Sobol indices for a reach of the Garonne river is presented.
The paper is organized as follows. Section 2 provides a brief overview of uncertainty in hydraulics. Section 3 presents the methods for dimension reduction, clustering and classification, and polynomial chaos for the mixture of experts surrogate. It also presents metrics to assess the validity of the surrogate and the formulation of Sobol indices. Results are presented in Sect. 4, illustrating the capability of the rMPCE to deal with both high-dimension and complex non-linear processes. Finally, concluding remarks, limitations, and perspectives are given in Sect. 5.
2 Uncertainty quantification for hydraulic modeling

The Garonne catchment
The study area extends over a 50 km reach of the Garonne river (southwest France) from Tonneins (upstream) to the confluence with the rivers Lot and La Réole (downstream) (see Fig. 1). It has a population of nearly 40,000 mainly concentrated in Tonneins and Marmande. This part of the valley is identified as an area at high risk of flooding (Lang and Coeur 2014 . Their characteristics are very different from one season to another, but the threats they represent remain very important. This part of the valley was equipped in the nineteenth century with infrastructure to protect the Garonne floodplain from flooding events. A system of longitudinal dykes and weirs was progressively constructed after the 1875 flood in order to protect floodplains and organize submersion and flood retention areas. Protections on the Garonne river form a system of successive storage areas for the floodplain beyond the dikes. This configuration is similar to the characteristic of other managed rivers such as the Rhone and the Loire. The QoI for the study is the water depth simulated over the river bed and the floodplain using the bidimension numerical model presented in Sect. 2.2. The uncertainties in the model parameters and forcing as well as in the model outputs are described in Sect. 2.3.

2D hydraulic modeling
The Shallow Water Equations (SWE) (de Saint-Venant 1871) are commonly used in environmental hydrodynamics modeling. They are derived from the Navier-Stokes equations (Sohr 2001) and based on the assumption that the horizontal length scale is significantly greater than the vertical scale, implying that vertical velocities are negligible, vertical pressure gradients are hydrostatic, and horizontal pressure gradients are due to displacement of the free surface. SWE express mass and momentum conservation averaged in the vertical dimension. The non-conservative form of the equations is written in terms of the water depth h and the horizontal components u x and u y of the velocity u ! in Cartesian coordinates (Hervouet 2007a): Momentum along x : Momentum along y : where grad ! and div are the gradient and divergence operators, with: To solve the system of SWE (1), initial conditions hðx; y; t ¼ 0Þ ¼ h 0 ðx; yÞ, u x ðx; y; t ¼ 0Þ ¼ u x;0 ðx; yÞ and u y ðx; y; t ¼ 0Þ ¼ u y;0 ðx; yÞ are provided along with boundary conditions (BC) at the surface, the bottom, and at upstream and downstream frontiers: hðx BC ; y BC ; tÞ ¼ h BC ðtÞ.
Due to the presence of non-linear terms in SWE, a closed-form solution of those equations is not available, except for very simplified cases. Therefore, they are discretized in space/time and their dynamic is numerically integrated using various schemes, e.g., method of characteristics (Chintu 1986), (discontinuous) Galerkin method (Eskilsson and Sherwin 2004), finite-element method (Hervouet 2007b), and finite-volume method (Anastasiou and Chan 1997), among others.
In this study, the Telemac-2D (T2D) 1 solver (Galland et al. 1991) based on a finite-element method is used (Hervouet 2007b). The equations are solved over a triangular mesh (see Fig. 1) featuring about 41,000 nodes, refined in the river bed and near the dykes. The discharge at Tonneins is imposed as the upstream boundary condition where the state-discharge rating curve at La Réole is imposed as the downstream boundary condition. A quasiunsteady state is considered, which refers to the convergence to a steady state. Indeed, the upstream discharge is set as a ramp starting from the initial condition value (1500 m 3 s À1 ) linearly increasing to a constant Q up (denoted Q for simplicity in the following). Each T2D transient simulation is integrated over 3 days (53 time steps of 5000 s) so that a steady flow associated to Q is prescribed over the entire area at the end of Day 3.
The Strickler friction coefficient K s is uniformly defined over four areas as displayed in Fig. 2. The friction coefficient values result from a calibration procedure over a set of non-over flowing events. These are set, respectively, to 45, 38, and 40 m 1=3 s À1 over upstream, middle and downstream parts of the river bed and 17 m 1=3 s À1 over the floodplain. More details on the Garonne river T2D model are given in Besnard and Goutal (2011).

Hydraulic uncertainty quantification
Typically, uncertainties are classified in two groups: epistemic uncertainty, resulting from incomplete knowledge of the correct settings of the model's parameters, and aleatory uncertainty, resulting from the incomplete knowledge of the true value of the physical system and usually linked to the aleatory nature of the physics. In this study, both epistemic and aleatory uncertainties are considered by investigating the effect of uncertainties in friction coefficients and in the upstream discharge forcing on water depth for the transient flow simulated with T2D.
Indeed, the small number of discharge and water depth measurements limits the spatial description and calibration of the friction in the river bed and the floodplain, leading to discontinuous values between friction areas. The K s coefficients setting is indeed prone to uncertainty related to the zoning assumption, the calibration procedure, and the set of calibration events. This uncertainty is more significant in the floodplain area where there is no observing station. The limited number of measurements also yields errors in upstream inflow to the river as it relies on the use of a rating curve, usually extrapolated for high flow, to translate the inflow from the measured water depth.
In the GSA sampling, the uncertainties in the friction coefficients and inflow are assumed to be independent. This assumption brings significant simplification with respect to reality where friction depends on water level. Yet it allows for a simplified description and calibration of friction coefficients, given the density of the observing network.
Classically, according expert knowledge, the friction coefficient is contained in an interval bounded by physical values depending on the roughness of soil material (Vazquez 2006;Goutal et al. 2018). Consequently, using the principle of maximum entropy (Shore and Johnson 1980), the distribution of the bounded Strickler friction coefficient is uniform. The boundaries of the uniform distribution are arbitrarily chosen AE5 from the calibrated value (Besnard and Goutal 2011) for the main channel roughness, as shown in Table 1. The Strickler friction coefficient of the floodplain is characterized by high uncertainty due to different land cover; therefore, the support of its distribution is wider and the boundaries have been chosen based on expert judgment. It should be noted that small Strickler's coefficient values are considered to account for the presence of vegetation or urban areas in the floodplain.
The upstream discharge is estimated using an extrapolation of discharge frequency curves at high probabilities (75 %) of occurrence of floods with a return period of two years. Confidence intervals on the extrapolated value can be derived. In that case, when the mean value (discharge of the two-year return period) and the standard deviation (extrapolated from the confidence intervals) are known, the maximum entropy distribution is Gaussian (Shore and Johnson 1980). The upstream discharge is, therefore, assumed to follow a Gaussian distribution centered on its biennial value at Tonneins (3300 m 3 s À1 ), with a standard deviation of 1100 m 3 s À1 . Moreover, to avoid unrealistic values, the PDF is truncated at 600 m 3 s À1 , corresponding to the annual mean discharge, and 6000 m 3 s À1 , corresponding to the vicennial flood at Tonneins. The characteristics of the uncertain model inputs distributions are summarized in Table 1. 3 Uncertainty propagation using reduced mixture of polynomial chaos expansions

Introduction to the rMPCE strategy
This section proposes a reduced Mixture of Polynomial Chaos Expansions (rMPCE). This advanced surrogate model strategy aims to predict a 2D output field subject to non-linearities with respect to sub-divided input space variables. This strategy features an output reduction stage and a local regression stage via clustering and classification. These stages are detailed in the following after a general presentation of the strategy. The direct model in denoted by M. It computes a p length real output y ¼ y 1 ; . . .; The learning set consists of n (input, output) samples, a.k.a., evaluations, snapshots, or observations, is denoted x ðiÞ ; y ðiÞ È É i2L , where L ¼ f1; . . .; ng is the set of indices of the n learning samples. The corresponding learning input matrix is denoted X with ½X ij ¼ x ðiÞ j and the learning output matrix is denoted Y with ½Y ij ¼ y ðiÞ j . Lastly, underline is reserved for random variables (e.g., u or U) while vectors and matrices are written in bold and in lower case (e.g., x) and upper case (e.g., X), respectively.
Transposed to the test case, these elements are defined as follows. The vector of upstream inflow and spatially defined friction coefficients x ¼ ðQ; K s Þ is denoted x when treated as a random variable. y ¼ h 1 ; . . .; h p À Á is the 2D water depth field at the T2D simulation time step of interest T, discretized over a mesh of size p and denoted y when treated as a random variable. The time step of interest corresponds to the flood' s rising part; it occurs 1 day, 2 h, 21 min and 20 s after the beginning of the studied flood. At this simulation time, the classical PCE leads to poor results (El Garroussi et al. 2020). Without loss of generality, the proposed strategy remains applicable for all time steps.
The rMPCE strategy is a two-stage process as illustrated in Fig. 3: 1. an offline learning stage that builds the model from a learning database, 2. an online prediction stage that evaluates the model to issue a prediction.  Moreover, the hyper-parameters of the surrogate model can be optimized in an outer loop around the learning stage in order to increase its accuracy measured on a validation database.
The learning stage developed in algorithm 1 features four main steps: 1. Reduction of the output variable dimension from p tõ p\p: the original space of dimension p is replaced with a latent space of dimensionp built from the learning output matrix Y. The learning output matrix Y 2 M n;p ðRÞ is then replaced with the reduced learning output matrixỸ 2 M n;p ðRÞ, which is computationally easier to handle. This reduction step is called encoding while the reverse is called decoding and maps from the latent space Rp to the original one R p . 2. Unsupervised clustering of the n learning output data into K groups, a.k.a., clusters: the reduced learning output matrixỸ is split into K local reduced learning output matrixỸ ðkÞ 2 M n k ;p ðRÞ; k 2 f1; . . .; Kg, where the n k observations inỸ ðkÞ share common patterns; L k & L is the sub-set of the learning indices of the samples belonging to the kth cluster, with [ K k¼1 L k ¼ L and L k \ L k 0 ¼ ; for any k 0 6 ¼ k.
3. Classification of the input space into K subspaces, based on the clustering results: • This step defines the boundaries of separation between the different classes within the input space. • This step provides a classifier taking a x as input and returning its degree of membership C k ðxÞ to the kth class, with C k ðxÞ ! 0 and P K k¼1 C k ðxÞ ¼ 1 by construction.
4. Construction of a 2D-functional output PCE surrogate for each cluster; e.g., for the kth cluster: • The dimension of the local output matrix Y ðkÞ ¼ y ðiÞ j i 2 L k 1 j p related to the kth cluster is reduced from p top and denoted g Y ðkÞ .
• The local surrogate model maps from the input space to the local latent space and requires a decoding step to go back to the original local output space. The prediction phase predicts the water depthŷ of a given input x. First, the degree of membership to the K classes is computed from the classifier: C 1 ðxÞ; . . .; C K ðxÞ. Then, the local PCE models are evaluated at x. Lastly, the global prediction in the latent space is a convex combination of the local ones: andŷ is expanded to the original output space, resulting in the predictionŷ. The current study is limited to hard classification, where a single class is attached to a given x. This implies that C k : R d ! f0; 1g instead of C k : R d ! ½0; 1. This results in the evaluation of a single local PCE; more precisely, the one indexed byk 2 fk : C k ðxÞ ¼ 1g.

Dimension reduction
In spite of recent advances that propose to estimate the PCE coefficients on a sparse grid (Eldred and Burkardt 2009) or with basis adaptive methods (Li and Ghanem 1998), the formulation of a surrogate model remains computationally expensive, especially when the dimension of the output is large. A common strategy applied here, is to build a surrogate model in a reduced output space, evaluating it for an input value, and then projecting its output value onto the original output space. In this study, two dimension reduction methods are investigated: PCA and AE. Both methods are applied on Y 2 M n;p ðRÞ, the matrix of the n evaluations of the p-length output y as illustrated in Fig. 3.
In this study, the output y is the water depth field discretized over the T2D unstructured mesh over the Garonne area. The output matrix Y is encoded onto a reduced latent space (see Fig. 4) as the matrix of the n evaluations of thep -length reduced outputỹ,Ỹ 2 M n;p ðRÞ, and is further used for the clustering stage. Moreover, any element of the latent space can be decoded onto the original output space. In particular, the initial matrix Y can be reconstructed, with some loss of information quantifying the performance of the reduction dimension technique.

Principal components analysis
PCA (Wold et al. 1987 The SVD of Y reads Y ¼ UDV > , where U is an n Â n orthogonal matrix, V is a p Â p orthogonal matrix, and D is a rectangular diagonal matrix with non-negative real numbers on the diagonal. Columns of UD are called principal components (PCs) and form an orthonormal basis in which the n samples y ð1Þ ; . . .; y ðnÞ are linearly uncorrelated. Then, the projection of the latter on thep p first PCs reads:Ỹ ¼ ½U :;1:p ½D 1:p;1:p 2 M n;p ðRÞ, thus reducing the output data dimension from p top. The column of V displays the corresponding weights associated to the PCs and any observationỹ in the latent space can be projected onto the original space: y ¼ V > Â Ã 1:p;1:pỹ . PCA allows summarizing data when the interesting patterns increase the variance of projections onto orthogonal components. But PCA also has limitations that are developed in Lever et al. (2017): the underlying structure of the data must be linear, patterns that are highly correlated may be unresolved because all modes are uncorrelated, and the goal is to maximize variance and not necessarily to find clusters.

Autoencoder
In order to deal with non-linear structure in the data matrix Y, the use of an AE (Hinton and Salakhutdinov 2006;van der Maaten et al. 2007) for dimension reduction was investigated. It relies on an unsupervised artificial neural network that encodes a variable of dimension p into a latent variable of dimensionp p and decodes this latent one to a recovered variable of dimension p, as close as possible to the original. The latent space is often called a bottleneck because of the particular shape of this neural network, illustrated in Fig. 4. In this paper, an AE with a symmetrical architecture (Nowlan and Hinton 1992) was used in order to reduce the number of parameters to be optimized in the network; it is based on encoder-decoder weight sharing. Steps needed for encoding and decoding are presented in Algorithms 3 and 4, respectively. An '-depth encoder maps y 2 R p onto the latent space Rp using ' successive encoding transformations: is a matrix of weight parameters, b l 2 R p l is a vector of bias parameters, and r l : R p l 7 !R p l is an activation function. The ' successive layers are of decreasing dimension: Then, the decoder maps the latent variable u ' 2 Rp onto the original space, using ' successive decoding transformations using the transposes of the encoder weight matrices as weight matrices for the decoder: with r 0 being the identity function. Thus, an autoencoder / ' ¼ / d;' / e;' is a sequence of 2' transformations, the first ' performing an encoding / e;' : R p 7 !Rp and the next ' a decoding / d;' : Rp7 !R p . The learning phase seeks the weights and biases minimizing the error k/ ' ðYÞ À Yk 2 2 while the use phase expands the dimension of a vectorỹ 2 Rp with the decoding function: / d;' ðỹÞ 2 R p . These weights and biases are usually initialized randomly and updated during training through the gradient backpropagation technique (Amari 1993).

Clustering and classification tools
Clustering is an unsupervised learning process that classifies data for which variables are observed via labels by using similarity measures. This approach is widely used for the purposes of data visualization, data compression, data denoising, or to better understand the correlations present in the data. In the present work, clustering methods are applied to the matrixỸ resulting from the output dimension reduction stage as shown in Fig. 3. It seeks to group the n dimension-reduced observationsỹ ð1Þ ; . . .;ỹ ðnÞ n o into K clusters and create the corresponding sub-sets of learning indices L 1 ; . . .; The kth cluster is associated with the label k, also known as the index or class. Then, these labels are mapped to the input space, here upstream forcing and bottom friction, to train a classifier x ! C 1 ðxÞ; . . .; C K ðxÞ ð Þ mapping from R d to ½0; 1 K to identify the boundaries between these clusters in the input space and to give the degree of membership of an input x to each of the corresponding classes. In the case of hard classification, only one class is associated to the input x and so the classifier maps from R d to f0; 1g K .

Clustering
Formally, clustering involves partitioning the set of observations L into K disjoint sets L 1 ; . . .; L K by returning labels indicating the index of the class of membership of each observation. Both k-means (Likas et al. 2003) and Gaussian mixture models (Mclachlan and Basford 1988) clustering algorithms are investigated in this paper. Both require prescribing the number of clusters K. The latter can either be prescribed manually based on the user's knowledge or estimated from a selection criteria such as the silhouette criterion (Rousseeuw 1987) that evaluates the separation distance between the resulting clusters. For a given observation indexed by i, belonging to the kth cluster, the silhouette criterion reads: where: • a k ðiÞ ¼ 1 jL k jÀ1 P j2L k ;j6 ¼i dði; jÞ is the average distance of the ith observation to all other observations in the kth cluster, with dði; jÞ ¼ỹ ðiÞ Àỹ ðjÞ 2 , • b k ðiÞ ¼ min l6 ¼k 1 jL l j P j2L l dði; jÞ is the smallest mean distance of the ith observation to all observations in any other cluster, of which i is not a member, • jSj is the cardinal number of the set S. Therefore, if i has been properly assigned, then the score s k is equal to 1. A score of 0 means that clusters are overlapping, and a score less than 0 means that i was assigned to the wrong cluster.
The k-means algorithm partitions the n observations into K clusters in which each observation belongs to the cluster with the nearest mean. It seeks to minimize the variance within the clusters: where l k ¼ jL k j À1 P i2L kỹ ðiÞ is the empirical mean ofỹ in cluster L k .
Given an initial set of K means l ð1Þ 1 ; :::; l ð1Þ K , the algorithm iterates the two following steps, presented at iteration t, until convergence: • Assignment step: assign each observation to the cluster with the nearest mean, i.e., 8k 2 f1; . . .; Kg: Because k-means struggles with clusters of varying density and with outliers, a clustering algorithm based on mixture of distributions was investigated here. The Gaussian Mixture Model (GMM) relies on the assumption that the empirical distribution of the n observed vectorsỹ ð1Þ ; . . .;ỹ ðnÞ is close to a mixture of K Gaussian distributions. Then, each observation is associated to the most likely Gaussian distributions, which then defines its cluster. The GMM is a mixture of K multivariate normal distributions. The kth distribution is characterized by its mean l k , covariance matrix R k , and weight x k . The PDF of this GMM reads: where p G ðỹ; l; RÞ ¼ 1 is the PDF of the Gaussian distribution with mean l and covariance matrix R. The GMM parameters x k ; l k ; R k f g 1 k K are estimated iteratively using an Expectation Maximization algorithm (Moon 1996;Bettebghor et al. 2011) until convergence of the likelihood. The expectation of the posterior probability k k of belonging to cluster L k can be expressed with Bayes' theorem: This is the E-step, where E stands for Expectation. Then, the mixture parameters l k and R k can be re-estimated by maximizing p G ðỹ; l k ; R k Þ: This is the M-step, where M stands for Maximization. The cluster of each observation i can be determined using Eq .4. Contrary to the k-means grouping, the GMM grouping can be either soft or hard. Soft grouping means that each observation i is assigned to each cluster in a weighted manner while hard grouping means that each observation i belongs to only one cluster. In this study, a hard splitting is considered: any pointỹ is assigned to cluster arg max k p G ðỹ; l k ; R k Þ.
The clustering assigns each of the n learning observations a label among f1; . . .; Kg. The classification uses these labels to draw the boundaries between classes in the input space.

Classification
Classification is a supervised learning process based on labels and derived from the clustering that groups observations into classes with respect to their labels, and identifies the boundaries between these classes.
The clustering has annotated each of the n learning observations with a label. According to these labels, the input variable, here, the liquid boundary condition and the friction of the river bottom x ðiÞ ¼ Q ðiÞ ; K ðiÞ s h i , is associated to kth cluster. The degree of membership of x ðiÞ to the kth cluster is written through the corresponding variable c i , such as c i ¼ k.
Here, a multi-class classification algorithm is considered: support vector machines (Cortes and Vapnik 1995). Support Vector Machines (SVM) aim at solving classification problems by finding good decision boundaries between two classes within the input space. For multiclass classification (K [ 2), the same principle is used. The multi-class problem is broken down to multiple binary classification cases called one-vs-one. SVM proceed to find the decision boundaries in two steps: • Mapping step: Input data are mapped to a new highdimension representation (target representation space) where the classification problem becomes simpler and where the decision boundary can be expressed as a hyperplane. • Maximizing the margin step: The separation hyperplane (decision boundary) is computed by maximizing the distance between the hyperplane and the closest data points from each class.
Because the mapping step is often computationally intractable, a ''kernel trick'' (Vapnik 1995;Scholkopf et al. 1999) is used. It is based on a kernel function k that maps any two input data x ðiÞ ; x ðjÞ È É to the distance between these data in the target representation space, completely bypassing the explicit computation of the new representation. The kernel trick is also used to develop non-linear generalization of the SVM. Let H be a k-kernels space. A general SVM is a discriminator of the form DðxÞ ¼ c i ðf ðxÞ þ bÞ where f 2 H and b 2 R are given by solving the general problem for a given C ! 0: . . .; ng; 0 f i ; 8i 2 f1; . . .; ng: where f i model the potential errors when the margin constraint is not verified. The decision functions of the following form are obtained: where A is the constraints set and a i are solutions of the following quadratic programming problem: 0 a i C; 8i 2 f1; . . .; ng; P n i¼1 a i c i ¼ 0: The main advantage of the SVM algorithm is its capability to deal with a wide variety of classification problems including high-dimension and non-linearly separable problems. One of its major drawbacks is that it requires many parameters to set correctly (under Scikit learn library (Pedregosa et al. 2011)) to attain good classification results.

Polynomial chaos expansions
The PCE surrogate model is built within each of the K classes in parallel (see Fig. 3). Let us consider the construction of a PCE within a single class and a computational model of interest M : D x & R d 7 !R, taking the vector x ¼ ðx 1 ; . . .; x d Þ 2 D x as input and returning y 2 R p as output: y :¼ MðxÞ. In the following, for the sake of simplicity, y is assumed to be a scalar (p ¼ 1). In the case of a vectorial response (p [ 1), the following derivations hold component-wise.
In uncertainty quantification, the deterministic input vector x is replaced by the associated random variable x ¼ Þand y ¼ MðxÞ is in turn a random variable. x is defined over the probability space ðD x ; F; PÞ and f x is its joint PDF. We seek to quantify the uncertainty in y due to uncertainty in x 1 ; . . .; x d . We assume that the random input variables are independent so as to comply with the assumption required for the polynomial chaos expansion theory. We also consider that the scalar output y is a second order random variable, i.e, E y 2 h i \ þ 1.
Under the previous assumptions, the random variable y can be expressed as a generalized polynomial chaos expansion (Xiu and Karniadakis 2002;Soize and Ghanem 2004): where w a ðxÞ ¼ Q d i¼1 w i;a i ðx i Þ is a tensor product of univariate orthonormal polynomials, i.e. E w i; c a is the deterministic coefficient associated with w a . a ¼ ða 1 ; :::; a d Þ is the multiindex vector with a i the degree of the univariate polynomial w i;a i and c a i ¼ y; w i a i ðxÞ Xiu and Karniadakis (2002) show the set of polynomials that provides an optimal basis for the different continuous probability distributions of the input variable x. It is derived from the family of hyper-geometric orthogonal polynomials known as the Askey scheme (Dongbin and Karniadakis 2003). The optimality of these basis selections derives from orthogonality with respect to weighting functions that correspond to the PDFs of the continuous distributions when placed in a standard form. For instance, when x i is a standard uniform (resp. standard normal) random variable, the corresponding basis comprises orthonormal Legendre (resp. Hermite) polynomials (Abramowitz et al. 1988).

Truncated polynomial chaos expansion
In practice, it is not tractable to use an infinite series expansion. An approximate representation is obtained with a truncation: with A 2 N d the truncation set of size m, i.e., c A ¼ c a ð Þ a2A 2 R m and A ðxÞ ¼ P a2N d nA c a w a ðxÞ the truncation-induced error. Blatman and Sudret (2011) introduced a hyperbolic truncation scheme that selects all polynomials satisfying the following criterion: with P being the highest total polynomial degree and 0\q 1 being the parameter determining the hyperbolic truncation surface. To further reduce the number of candidate polynomials, one can additionally apply a low-rank truncation scheme that reads (Sudret 2015): where kak 0 is the rank of the multivariate polynomial w a , defined as the total number of non-zero components a i ; i ¼ 1; :::; d. In this study, the prescribed rank r is chosen as a small integer value, e.g., r ¼ 2; 3 (Mai et al. 2016) and the polynomial degree P is varied from 2 to 9, and the value retained is the one that minimizes the prediction error.

Estimation of coefficients
The computation of the coefficients c a in Eq. 9 can be conducted by means of intrusive (i.e, Galerkin scheme) or non-intrusive approaches (e.g., stochastic collocation, projection, regression methods) (Blatman et al. 2007). In this paper, we consider a standard regression method based on the minimization of a mean squared learning error (-Baudin et al. 2017). In practice, the coefficients are obtained by minimizing an empirical mean over a learning database: where x ðiÞ ; i 2 L È É is a Design Of Experiment (DOE) obtained with a random sampling of the input random vector. For that purpose, the computational model M is integrated for each point of the DOE, yielding the learning output matrix Y. Equation 10 basically represents the problem of estimating the parameters of a linear regression model, for which the least squares solution readŝ is the information matrix containing the evaluation of the polynomial basis functions over the DOE. Hence, the approximated output variableŷ can be expressed as follows: At the prediction phase, only the PCE related to the class to which the new observation belongs is evaluated (hard evaluation).

Surrogate model validation metrics
In the present study, two standard metrics are used to measure the quality of the rMPCE surrogate model at T: the Q 2 predictive coefficient and the Root Mean Squared Error (RMSE). The validation is carried out over an (input, output) validation database D v of size n v .

Predictive coefficient
At the kth mesh node, the Q 2 predictive coefficient is defined as: where MSE k ðD v Þ ¼ n v À1 P n v i¼1 y ðnþiÞ k Àŷ ðnþiÞ k 2 and MSE k ðD v ; meanÞ ¼ n À1 v P n v i¼1 y ðnþiÞ k À y k 2 is the MSE of the averaging model returning the mean of the learning outputs whatever the input parameter value. The global counterpart of MSEðD v ; meanÞ is computed spatially by averaging over the p elements of the output vector: Thus, the global counterpart of Q 2 is: The predictive coefficient measures the performance of the surrogate model with respect to the data average. When Q 2 is lower than (resp. equal to) zero, the surrogate is worse than (resp. equal to) the learning output values average. When Q 2 is equal to one, the surrogate interpolates the validation database. In practice, the surrogate is deemed appropriate when Q 2 is greater than 0.8. The predictive coefficient is also found under the name of Nash-Sutcliffe model efficiency coefficient in the hydrological literature, where it assesses the predictive capacity of the simulated discharge over a time window with respect to observed discharges (Nash and Sutcliffe 1970).

Root Mean Squared Error
The RMSE is used to measure the accuracy of the model and should be equal to 0 when the model is perfect. At the kth given mesh node, it is defined as the square root of the mean squared errors (MSE), measuring the squared distance between the surrogate model and the reference model: Their global counterpart are:

Sensitivity analysis
Sensitivity analysis aims to investigate how the different uncertain input variables x 1 ; . . .; x d influence the output variable y ¼ MðxÞ over the whole uncertain input space. M either stands for the direct solver or for its surrogate. The overall objective is to identify which input parameters contribute the most to the uncertainty in the output and to order them accordingly. For the sake of simplicity, we focus on a mono-dimensional output variable y. The model output uncertainty can be represented by its variance V½y to be explained on the basis of the uncertain input variables and their interactions. This is the purpose of the Sobol methodology (Sobol 1993;Saltelli 2010;Iooss and Lemaître 2015;Razavi et al. 2021), valid when x 1 ; . . .; x d are independent and when y is a second-order random variable, i.e, E y 2 h i \1. This technique decomposes the total output variance V y h i into 2 d À 1 elementary contributions: where: • I d ¼ f1; . . .; dg; is the contribution of the x i in interaction with x j ; • and so on.
In practice, interest is focused on standardized versions of these contributions: X is the Sobol index related to the interaction between the uncertain input variables x i ; i 2 u. S u is the part of V y h i explained by this interaction. All these indices add up to 1 and, thus, represent proportions of output variance. Most of the time, Sobol study is conducted on: • the first-order indices, S 1 ; . . .; S d , where S i represents the part of V½y explained by x i only; and • the total-order indices, S T 1 ; . . .; When the difference between S i and S T i is significant, this means that there are interactions between x i and other uncertain input variables explaining V½y. In this case, it is common to look at the value of the second-order indices S i;1 ; . . .; S i;d and so on. Conversely, S T i % S i leads to the conclusion that there is no interaction between x i and another variable explaining V½y. Consistently, P i S i ¼ 1 if there is no interaction between the input parameters.

Strategy and experimental settings
The rMPCE strategy results at T are compared to those of a classical PCE strategy. Different choices for dimension reduction, clustering, and regression are investigated. For this purpose, two databases are generated in this study with an optimized Latin Hypercube Sampling (LHS) (Damblin et al. 2013) for the uncertain input variables whom PDFs are described in Table 1: • a learning database of 1000 T2D evaluations to build and fit the surrogate model; and • a validation database of 500 T2D evaluations to evaluate the accuracy of the surrogate model.

Computational environment
CERFACS's cluster, Nemo, has been used to run T2D simulations. The Nemo cluster includes 6912 cores distributed in 288 compute nodes. The ECU power peak is 277 Tflop/s. The computational cost of T2D solver is reduced thanks to the parallel computing (single simulation lasts 6 min using 24 processors instead of 20 min using one processor). GSA based on a large set of T2D simulations is too costly. Hence the need for surrogate model formulation. The rMPCE surrogate model proposed in this study is based on algorithms from different Python libraries. The first step uses AE from Keras Tensorflow (Géron 2017) with a graphics processing unit (GPU) support Python package to reduce the dimension of the output space. The second step of this algorithm involves clustering and classifying data using a GMM and SVM algorithms from the Scikit-Learn library (Pedregosa et al. 2011). In the final step, the algorithm constructs a local regression model within the cluster; for this purpose, PCE of the Open-TURNS library (Baudin et al. 2017) is used.
The meta-model learning stage (see Algorithm 1) is moderately costly: the tuning of the AE parameters takes about 3 h and the construction of the PCE takes about 15 min. The computational cost of the prediction stage is then drastically reduced, e.g., predicting 500 simulations takes 470 s.

Output dimension reduction
Dimension reduction results for PCA and AE are presented in Fig. 5. The size of the latent spacep is plotted along the x-axis, the left y-axis represents the RMSE (quadratic error between initial and reconstructed water level field) in meters for PCA (solid blue line) and AE (dotted blue line), and the right y-axis represents the cumulated explained variance for the PCA. Different neural network architectures were tested in order to minimize the RMSE metric. The resulting neural network is compiled with mean squared error loss and Adam optimizer (Zhang 2018) with 0.001 learning rate and the default Keras parameters. The number of training epochs is set to 200 while the batch size for the training cycle is set to 50. The size of the input is set to 41,416 neurons corresponding to the number of features in the database.
For PCA, the RMSE decreases exponentially from 9 to 3.82 centimeters as the number of principal components in the latent space increases from 1 to 50. For 26 components, 98% of the variance of the water depth is explained and the RMSE is about 4 centimeters. For a small number of components, AE leads to a larger RMSE than PCA: 27 centimeters against 9 centimeters for a single component. Beyond 24 components, AE leads to a smaller RMSE than PCA: 1.27 centimeter against 3.82 centimeters for 50 components. A latent space spanned over 37 components offers a good compromise between accuracy and computational cost for both methods. Despite the fact that AE is relatively expensive compared to PCA (2 h against 3 min), it allows to account for non-linearities in areas with strong gradient bathymetry, mainly in ditches and downstream of dikes. Indeed, the maximum absolute error for water depth reconstructed from the PCA displayed in Fig. 6 reaches 3 meters in a mesh node located in a ditch for a selected simulation, while the maximum absolute error for water depth reconstructed from AE remains smaller than 1 centimeter. Therefore, in the following, dimension reduction is achieved using the more accurate AE technique. Figure 7 displays the silhouette criterion defined in Eq. 2 for both k-means (top panels) and GMM (bottom panels) clustering methods, setting the number of clusters to K ¼ 2; 3; 4 (from left to right). The silhouette criterion s k ðiÞ is plotted along the x-axis for each observation i. The observation labels are indicated along the y-axis and arranged by the color-coded cluster number. The red vertical line indicates the average silhouette criterion computed among all observations and all clusters. This figure displays the quality of the clustering as well as the size of the resulting clusters. When K ¼ 2 and K ¼ 4, the size of the clusters are heterogeneous with silhouette values s k ðiÞ smaller than the mean value. K ¼ 3 provides homogeneous clusters with satisfying silhouette values for all clusters. In the following, the three classes resulting from the GMM classification are kept.

Clustering and classification
A hydraulic analysis of the clusters shows that the first cluster gathers medium-flow simulations where the flow submerges the dikes and barely propagates in the floodplain. The second cluster characterizes high-flow simulations where the flow significantly propagates in the floodplain and the third cluster characterizes low flow simulations where the flow is confined in the river bed.
The first (top panels) and second AE modes (bottom panels) for each cluster (with K ¼ 3) are shown in Fig. 8. In cluster 1, the first mode represents the mean flow dynamics while the second mode represents the flow obstacles. In cluster 2, the first mode corresponds to the maximum extent of the water while the second mode highlights the influence areas of upstream and downstream boundary conditions. In cluster 3, the first mode could be interpreted as the maximum flow extent and the second mode as a versus upstream-downstream flow. Figure 9 displays the predictive coefficient when the flood occurs, computed between the validation database and a classical PCE prediction on the left panel and between the validation database and the rMPCE prediction on the right panel. The areas where Q 2 is close to 0 are indicated in yellow and it clearly appears that rMPCE provides a far more predictive surrogate than classical PCE.

Regression
The classical PCE poorly predicts 6625 nodes (Q 2 \0:8) out of the 41,416 mesh nodes, mostly located in the floodplain where the response in water depth to change in friction and in inflow is non-linear ( Fig. 9 left panel). The rMPCE leads to a significant improvement for 90% of Fig. 5 Evolution of the RMSE computed between the real water depth (learning database) and the one reconstructed with the PCA inverse method in solid blue line and the AE decoder in dashed blue line, and of the reconstructed output variance for the PCA in solid red line, according to the latent space dimensionp Fig. 6 Spatialized maximum absolute error computed between a simulated water depth (one simulation from the learning database) and its reconstruction using PCA inverse method with 37 principal components. In zoom, the bathymetry profile along the horizontal section including the point with the maximum reconstruction error these poorly predicted nodes (564 nodes remain with Q 2 \0:8) as illustrated in Fig. 9 (right panel). Given these two maps, the contribution of the rMPCE strategy is significant for a good prediction of the water height in the floodplain where human and economic stakes are predominant. A zoom on the diagnostic of the rMPCE strategy for poorly predicted nodes (Q 2 \0:8) raises the contribution of the loop on the polynomial degree P to quality prediction improvement. The Q 2 resulting from setting the same P for the different classes in rMPCE (uniform rMPCE) is plotted as a dotted blue line in Fig. 10. P equal to 4 for the local PCEs of the three classes returns a value of Q 2 equal to 0.64, which is physically unsatisfactory. P greater than 4 leads to an over-fitting of the model to the learning data and a lower value leads to an underfitting.
The Q 2 resulting from the polynomial degree optimization loop (varying the polynomial degree P between 2 and 9) for PCE in each of the three classes is plotted as a dashed red line. The first class, mostly defined by medium flows, requires a P equal to 5 in order to approximate properly the water depth while the second class, characterized by high flows, requires a P equal to 4 and the third class, corresponding to low flows, requires a P equal to 3. This suggests that the physics in the first class is complex and requires to increase P, whereas the physics in the third class is rather simple as the optimal P for this class is equal to 3. Thus, the PCE polynomial degree optimization loop allows obtaining a good approximation of the modes representing the water depth following dynamics within each class.

Sensitivity analysis
The variance-based GSA in this study is based on Saltelli's method for the estimation of Sobol indices using the rMPCE surrogate model. The main goal of GSA is to rank the uncertain parameters according to their influence on the variance of the QoI, here, the water depth 2D field. Figure 11 displays the first-order (left panels) and total order (right panels) Sobol indices for the four Strickler friction coefficients and discharge (from top to bottom) at time T. Analysis of the first order Sobol indices reveals the large influence of the discharge as this uncertain variable explains about 80% of the water depth variance on the overall domain. The Strickler friction coefficient associated to the floodplain area influences by 9% the water depth variance upstream and in some dyked areas. The influence of the Strickler coefficients associated with the river bed remains weak or slightly significant in a few places; for example, K s;4 influences the water depth variance by 82% locally in a dyked zone downstream of the river. Fig. 9 Spatialized predictive coefficient computed between the validation database and the surrogate prediction at T: classical PCE (left) and rMPCE (right) Fig. 10 Evolution of the predictive coefficient Q 2 of rMPCE (solid red line) in which the polynomial degree has been optimized within each class and of uniform rMPCE (solid blue line with cross marker), resulting from setting the same polynomial degree for the different classes Fig. 11 Sobol indices of the hydraulic input variables estimated using Saltelli's method based on rMPCE for the simulated water depth at time T = 95,000 s. First-order indices are plotted on the left panels and total order on the right panels for K s;1 (floodplain), K s;2 (upstream river bed), K s;3 (middle river bed), K s;4 (downstream river bed), and Q (upstream forcing) from top to bottom The analysis of the total Sobol indices indicates that while the friction coefficients have a low first order Sobol index, they are not negligible as they have a significant influence through their interactions with other variables. Yet, the discharge remains by far the most influencing variable when it interacts with the other variables as shown in the right-bottom plot. It should be noted that the GSA results depend on the hypothesis on the input random variables distributions. For instance, the significant influence of the floodplain Strickler friction coefficient compared to that of the river bed coefficients may be due to the large uncertainty translated by the large range of K s;1 's uniform distribution.
5 Conclusions, limitations, and future research

Conclusions
In this paper, an rMPCE surrogate model is used to conduct a GSA in order to rank the sources of uncertainty with a variance-based sensitivity analysis in the presence of nonlinearities and at a parsimonious computational cost. The rMPCE strategy is based on a mixture of a polynomial chaos expansions implemented in a reduced output space and into clusters where non-linearities between input and output remain small. It is used to approximate the 2D water depth simulated using the T2D numerical solver. The uncertain input space contains five scalars and the uncertain output space is a 2D discretized field of large dimension (about 41,000 mesh nodes). This strategy is illustrated when the flood front enters the floodplain, causing nonlinearities between inflow, friction and the water field, especially in regions of strong bathymetry gradient. The first step of the rMPCE strategy involves compressing the water depth data. To this end, the PCA and AE methods were compared. PCA is a simple linear transformation on the input space to directions of maximum variation while AE is an advanced technique that minimizes the reconstruction loss. The AE technique yielded more accurate results as it was able to deal with non-linearities in the output field.
The second step of the rMPCE strategy involves grouping the reduced data with similar patterns into classes. After comparing the silhouette coefficient derived from the k-means and GMM methods, three classes were considered based on the GMM, leading to three different hydraulic behaviors. The third step consists of defining the boundaries between these classes within the input space using the SVM algorithm. It appears that the boundaries were mostly driven by the discharge variable.
The last step of the rMPCE strategy is to construct a local optimized PCE within each class. It was shown that the resulting surrogate model simulates properly the water depth over the study area and improves the prediction by 90% compared to the one given by a classical PCE. Indeed, PCE was successful in predicting water depth for over 83% of the grid points, mostly in the river bed. However, it fails to predict water depth in the floodplains where non-linearities occur. In these regions, rMPCE was able to deal with non-linearities and provide good prediction for 98% of the grid points.
Sobol indices were then estimated using the rMPCE surrogate model. It was shown that the water depth over the considered study area is predominantly controlled by the upstream discharge except for the left bank side of the upstream which is influenced by the Strickler friction coefficient of the floodplain. The total Sobol indices of the three Strickler friction coefficients related to the river bed indicate that despite the fact that those variables have a low first-order Sobol index in all domain, they are not negligible as they influence the water depth through interactions with the other variables. It has also been emphasized that those results depend on the description of the input variables PDF.

Limitations
In practice, tuning the AE hyper-parameters, such as the number of layers and the number of neurons per layer, remains difficult (van der Maaten et al. 2007). One way to overcome this limitation is to consider an existing architecture that was proven successful for a similar problem, or training the AE directly from the PCA response given that the AE may be considered as a non-linear extension of PCA, or using pre-training methods allowing for a layerby-layer learning (Makhzani and Frey 2015).
Due to time constraints, the model has not been tested for the case where all time steps are taken into account in one batch. Eventually, this could reduce the non-linearities present, in particular for the dimension reduction step.
The assumption for the description of the PDF for the Strickler coefficients could be revisited. An ensemble of coupled sediment-hydrology simulations could be generated in order to investigate how the topography evolves with the flow and consequently the friction evolves.
The assumption of independence of the input variables, Strickler's friction coefficient and upstream discharge, can be reviewed. In this sense, a sensitivity analysis could be conducted by considering the Shapely indices (Iooss and Prieur 2017).

Future research
As a perspective, first, the proposed surrogate modeling strategy should be applied to all time steps of the hydraulic simulation and to the computation of the time-varying Sobol indices. Also, numerical improvement could be reached with the analytical computation of Sobol indices from the local polynomial coefficients instead of their stochastic estimation with the Saltelli method with the rMPCE surrogate as implemented here. Additionally, the mixture strategy could also be revisited with kernel basedclustering methods that could take into account the nonlinearities, an adaptive re-sampling in clusters with a small predictive coefficient, and a weighted sum of the predictions from the local models using frequentist model averaging or Bayesian model averaging. A local mesh refinement in areas where the predictive coefficient of rMPCE remains small could be investigated. This would lead to further improvement relying on multi-fidelity approaches.
Another perspective would be to improve the rMPCE to simulate the hydraulic state on another time window than the one used for training in order to better meet the needs of data assimilation, typically when we go from one assimilation cycle to another. In this sense, a possible approach would be to combine rMPCE with NARX (Mai et al. 2016) to simulate the dynamics from one time step to another.
A major perspective for this work is to extend the uncertain input space. To begin with, the input space could include time-varying upstream forcing in order to simulate realistic flood events. It could also include a spatially refined friction field, potentially resulting from calibration with a densified, remotely sensed observation network. In both cases, the dimension of the input space should also be reduced, for instance, using the dimension reduction techniques applied here for the output space dimension reduction.
Finally, the resulting surrogate model can be used in the context of data assimilation. Indeed, the computation of the Sobol indices allows the identification of variables that should be included in the control vector. Then, the surrogate model could be used in place of the direct numerical solver for a low-cost stochastic estimation of the background covariance matrix in ensemble-based data assimilation algorithms. The assimilation of in-situ and remotesensing water level data with a parsimonious ensemblebased algorithm paves the way for the improvement of forecasted water depth and discharge in an operational framework.
Acknowledgements Funding for this work was provided by CER-FACS and Region Occitanie. This work was partly supported by the French national program LEFE/INSU. The authors gratefully thank the Electricité de France (EDF) for providing the Telemac 2D model for the Garonne river. The authors would like also thank R. Lebrun, C. Lapeyre, and S. Boyaval for their expertise and for their fruitful discussions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.