Comparison of Polynomial Chaos and Gaussian Process surrogates for uncertainty quantification and correlation estimation of spatially distributed open-channel steady flows

TL;DR

This paper compares Polynomial Chaos and Gaussian Process surrogates for efficient uncertainty quantification in steady open-channel flows, demonstrating their accuracy and computational advantages over traditional Monte Carlo methods.

Contribution

It introduces and evaluates surrogate models based on Polynomial Chaos and Gaussian Processes for rapid and accurate uncertainty estimation in river flow simulations.

Findings

Both surrogates match Monte Carlo accuracy with fewer simulations.

A small number of model evaluations (10-100) suffices for reliable covariance and PDF estimates.

Surrogates are effective even in complex flow regions with heterogeneous bathymetry.

Abstract

Data assimilation is widely used to improve flood forecasting capability, especially through parameter inference requiring statistical information on the uncertain input parameters (upstream discharge, friction coefficient) as well as on the variability of the water level and its sensitivity with respect to the inputs. For particle filter or ensemble Kalman filter, stochastically estimating probability density function and covariance matrices from a Monte Carlo random sampling requires a large ensemble of model evaluations, limiting their use in real-time application. To tackle this issue, fast surrogate models based on Polynomial Chaos and Gaussian Process can be used to represent the spatially distributed water level in place of solving the shallow water equations. This study investigates the use of these surrogates to estimate probability density functions and covariance matrices at a reduced computational cost and without the loss of accuracy, in the perspective of ensemble-based data assimilation. This study focuses on 1-D steady state flow simulated with MASCARET over the Garonne River (South-West France). Results show that both surrogates feature similar performance to the Monte-Carlo random sampling, but for a much smaller computational budget; a few MASCARET simulations (on the order of 10-100) are sufficient to accurately retrieve covariance matrices and probability density functions all along the river, even where the flow dynamic is more complex due to heterogeneous bathymetry. This paves the way for the design of surrogate strategies suitable for representing unsteady open-channel flows in data assimilation.

Tables

Table 1. Table 1: Q 2 subscript 𝑄 2 Q_{2} error for pGP and PC surrogates computed with respect to the MC reference. The error corresponds to the average over the M 𝑀 M stations with increasing number of snapshots N 𝑁 N from 49 to 256.

$N$	pGP	PC
49	0.99965	0.99983
121	0.99514	0.99993
256	0.99143	0.99962

Table 2. Table 2: Two-sample Kolmogorov-Smirnov statistical test for pGP and PC surrogates with respect to the MC reference at Marmande with increasing number of snapshots N 𝑁 N (49, 121, 256). The null hypothesis is rejected if D > 6 , 082 ⋅ 10 − 3 𝐷 6.082e-3 D>\numprint{6.082e-3} .

Surrogate	Snapshots	Statistics D	p-value
	49	$7,95\cdot {10}^{-3}$	0.004
pGP	121	$3,97\cdot {10}^{-3}$	0.409
	256	$3,02\cdot {10}^{-3}$	0.751
	49	$7,15\cdot {10}^{-3}$	0.012
PC	121	$4,95\cdot {10}^{-3}$	0.172
	256	$4,93\cdot {10}^{-3}$	0.175

Equations

$$\mathbf{x}\in {\mathbb{R}}^d\to \mathbf{h}=\mathcal{M}(\mathbf{x})\in {\mathbb{R}}^M.$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle{\partial_{t}A(h)+\partial_{a}Q=0}\\ \displaystyle{\partial_{t}Q+\partial_{a}\left(\frac{Q(h)^{2}}{A(h)}\right)+g\,A(h)\,\partial_{a}h-g\,A(h)\,(S_{0}-S_{f})=0}\\ \end{array}\right.$$

$$S_f=\frac{Q^2}{K_s^{2}A(h{)}^{2}R(h{)}^{4/3}},$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle{\partial_{a}Q=0}\\ \\ \displaystyle{\partial_{a}h=\frac{(S_{0}-S_{f})}{1-Fr^{2}}}\\ \end{array}\right.$$

$$F{r}^2=\frac{Q^2}{g{A}^3}\frac{\partial A}{\partial h}.$$

$${\hat{h}}_a\left(\mathbf{x}\right)=\sum _{i=0}^r{\gamma }_{a,i}{\Psi }_i\left(\mathbf{x}\right),$$

$${\int }_Z{\Psi }_i(\mathbf{\zeta }){\Psi }_j(\mathbf{\zeta })\rho (\mathbf{\zeta })\mathrm{d}\mathbf{\zeta }={\delta }_{ij},$$

$$r_{\text{pc}}=\frac{(d+P)!}{d!P!},$$

$${\hat{h}}_{\text{pc},a}(\mathbf{x}):={\mathcal{M}}_{\text{pc},a}(\mathbf{\zeta })=\sum _{\genfrac{}{}{0ex}{}{i\in {\mathbb{N}}^d}{|i|\le P}}{\gamma }_{a,i}{\Psi }_i\left(\mathbf{\zeta }\right).$$

$${\gamma }_{a,i}=<h_a,{\Psi }_i>\cong \sum _{k=1}^N{h}_a^{(k)}{\Psi }_i({\mathbf{\zeta }}^{(k)})w^{(k)},$$

$$\mathbf{Y}=\mathbf{U}\mathbf{\Lambda }{\mathbf{V}}^{\mathrm{T}}=\sum _{k=1}^{r_p}{\lambda }_k{\mathbf{u}}_k{\mathbf{v}}_k^{\mathrm{T}},$$

$$h_a({\mathbf{x}}^{(k)})=(\mathbf{U}\mathbf{\Lambda }{\mathbf{V}}^{\mathrm{T}}{)}_{ak}=U_{a:}(\mathbf{\Lambda }{\mathbf{V}}^T{)}_{:k}=\sum _{i=1}^{r_p}{\gamma }_{a,i}{\Psi }_i({\mathbf{x}}^{(k)}),$$

$${\hat{h}}_{\text{pgp},a}({\mathbf{x}}^{\ast })=\sum _{i=1}^{r_p}{\gamma }_{a,i}{\Psi }_{\text{gp},i}({\mathbf{x}}^{\ast }),$$

$${\tilde{\Psi }}_i(\mathbf{x})\sim \text{GP}\left(0,{\sigma }_i^2{\pi }_i(\mathbf{x},{\mathbf{x}}^{\prime })\right),$$

$${\pi }_i(\mathbf{x},{\mathbf{x}}^{\prime })=\exp \left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}^{\prime }{\Vert }^2}{2{\ell }_i^2}\right),$$

$${\Psi }_{\text{gp},i}({\mathbf{x}}^{\ast })=\sum _{k=1}^N{\beta }_{k,i}{\pi }_i\left({\mathbf{x}}^{\ast },{\mathbf{x}}^{(k)}\right),$$

$${\mu }_{h_a}$$

$${\sigma }_{h_a}$$

$$\mathbf{C}=\frac{1}{N_{\text{ref}}-1}\sum _{k=1}^{N_{\text{ref}}}\left({\mathbf{h}}^{(k)}-\overline{\mathbf{h}}\right){\left({\mathbf{h}}^{(k)}-\overline{\mathbf{h}}\right)}^T,$$

$${\mu }_{h_{\text{pc},a}}$$

$${\sigma }_{h_{\text{pc},a}}$$

$${\mathbf{C}}_{\text{pc}}={\left(\text{cov}\left(h_{pc,a_m},h_{pc,a_n}\right)\right)}_{1\le m,n\le M},\text{cov}\left(h_{pc,a_m},h_{pc,a_n}\right)=\sum _{\genfrac{}{}{0ex}{}{i\in {\mathbb{N}}^d}{|i|<P,i\ne 0}}{\gamma }_{a_m,i}{\gamma }_{a_n,i},$$

$$\mathbb{V}[{\hat{h}}_a]=\sum _{i=1}^d{\mathbb{V}}_i({\hat{h}}_a)+\sum _{j=i+1}^d{\mathbb{V}}_{ij}({\hat{h}}_a)+\cdots +{\mathbb{V}}_{1,2,...,d}({\hat{h}}_a),$$

$$S_i^a=\frac{{\mathbb{V}}_i({\hat{h}}_a)}{\mathbb{V}({\hat{h}}_a)},S_{ij}^a=\frac{{\mathbb{V}}_{ij}({\hat{h}}_a)-{\mathbb{V}}_i({\hat{h}}_a)-{\mathbb{V}}_j({\hat{h}}_a)}{\mathbb{V}({\hat{h}}_a)}.$$

$$S_{T_i}^a=\sum _{\genfrac{}{}{0ex}{}{I\subset \{1,\dots ,d\}}{I\ni i}}{S}_I^a.$$

$$S_{\text{pc},i}^a=\frac{1}{{\sigma }_{h_{\text{pc},a}}^2}\sum _{\genfrac{}{}{0ex}{}{j\in {\mathbb{N}}^d\text{s.t.}|j|<P,}{j_i>0\text{and}{j}_{k\ne i}=0}}{\gamma }_{a,j}^2,$$

$$Q_{2,a}=1-\frac{\sum _{k=1}^{N_{\text{ref}}}{\left(h_a^{(k)}-{\hat{h}}_a^{(k)}\right)}^2}{\sum _{k=1}^{N_{\text{ref}}}{\left(h_a^{(k)}-{\overline{h}}_a\right)}^2}.$$

$$D=\underset{x}{\sup }|F_n(x)-G_m(x)|.$$

$$D>c(\alpha )\sqrt{\frac{n+m}{nm}},$$

$$\text{RMSE}=\sqrt{\frac{1}{M}\sum _{a=1}^M({\square }_a-{\square }_{\text{mc},a}{)}^2.}$$

Full text

††journal: Stochastic Environmental Research and Risk Assessment

Nomenclature

• $\boldsymbol{\zeta}$

Standardized random variable

• $\delta$

Kronecker delta function

• $\ell$

Correlation length scale

• $\gamma$

Basis coefficient

• $\lambda$

Singular value

• $\mathbb{E}$

Expectation operator

• $\mathbb{V}$

Variance

• $\mathbf{\Lambda}$

Rectangular singular value matrix

• $\mathbf{\Pi}$

Correlation matrix evaluated for the database $\mathcal{D}_{N}$

• $\mathbf{C}$

Snapshot covariance matrix

• $\mathbf{h}$

Output random vector of size $M$

• $\mathbf{U}$

Orthogonal square left singular matrix

• $\mathbf{V}$

Orthogonal square right singular matrix

• $\mathbf{x}$

Input random vector of size $d$

• $\mathbf{x}^{*}$

Input random vector not included in the training set $\mathcal{X}$

• $\mathbf{Y}$

Centred snapshot matrix

• $\mathcal{D}_{N}$

MASCARET simulation database of size $N$

• $\mathcal{M}$

MASCARET forward model operator

• $\mathcal{N}$

Normal distribution characterized by mean and STD

• $\mathcal{U}$

Uniform distribution characterized by minimum and maximum

• $\mathcal{X}$

Input training set of size $N\times d$

• $\mathcal{Y}$

Output training set of size $N\times M$

• $\mu$

Mean value

• $\omega$

Gaussian quadrature weight

• $\pi$

Correlation function

• $\Psi$

Basis function

• $\rho$

Joint probability density function

• $\sigma$

STD

• $\tau$

Nugget effect

• $\widehat{h}$

Estimated water level

• $a\in[a_{\text{in}};a_{\text{out}}]$

Curvilinear abscissa $\text{\,}\mathrm{k}\mathrm{m}$

• $A$

Hydraulic section $\text{\,}\mathrm{m}^{2}$

• $c(\alpha)$

$\alpha$-level tabulated value associated with $D$

• $D$

Kolmogorov-Smirnov statistic

• $d$

Uncertain space size

• $h$

Water level part of the hydraulic state $(h,Q)$ $\text{\,}\mathrm{m}$

• $i$

Surrogate decomposition index

• $k$

Snapshot index

• $K_{s}$

Strickler friction coefficient $\text{\,}\mathrm{m}^{1/3}\,\mathrm{s}^{-1}$

• $M=\numprint{14}$

Number of observation stations

• $N$

Training set size

• $N_{\text{ref}}=\numprint{100,000}$

Validation set size

• $P$

Wet perimeter $\text{\,}\mathrm{m}$

• $Q$

Discharge part of the hydraulic state $(h,Q)$ $\text{\,}\mathrm{m}^{3}\,\mathrm{s}^{-1}$

• $Q_{2}$

Predictive squared correlation coefficient

• $R$

Hydraulic radius $\text{\,}\mathrm{m}$

• $r$

Number of terms in surrogate model

• $RC$

Upstream water level-discharge local rating curve

• $S_{0}$

Channel slope $\text{\,}\mathrm{m}\,\mathrm{k}\mathrm{m}^{-1}$

• $S_{f}$

Friction slope

• $S_{i}$

First order Sobol’ index

• $S_{T}$

Total order Sobol’ index

• $W$

River width $\text{\,}\mathrm{m}$

• CDF

Cumulative Distribution Function

• DA

Data Assimilation

• EnKF

Ensemble Kalman Filter

• GP

Gaussian Process

• MC

Monte Carlo

• PC

Polynomial Chaos

• PDF

Probability Density Function

• PF

Particle Filter

• PGP

POD-based Gaussian Process

• POD

Proper Orthogonal Decomposition

• RMSE

Root Mean Square Error

• SA

Sensitivity Analysis

• STD

STandard Deviation

• SVD

Singular Value Decomposition

• SWE

Shallow Water Equation

• UQ

Uncertainty Quantification

1 Introduction

The predictive skills of hydraulic models have greatly increased with advances in free surface flow numerical modelling and computational resources. Real-time flood forecasting relies on the use of sparse in situ observations as well as imperfect hydrology or hydraulic models usually solving the 1-D Shallow Water Equations (SWE). Assessing the predictive capabilities of these hyperbolic partial differential equations remains an important challenge as public safety and water resource management are at stake (Weerts et al. 2011). SWE solve for spatially varying water level and river discharge (referred to as the river state) using physical parameters (e.g. friction coefficients), initial conditions and boundary conditions described as a hydrograph or limnigraph. These input data are subject to epistemic uncertainties due to an imperfect knowledge of the river properties as well as to aleatory uncertainties related to environmental and meteorological intrinsic hazards. Both types of errors translate into uncertainties in the simulated river state, thus preventing the hydraulic model from being effective in forecast mode. In practice, these uncertainties can be reduced when complementary data become available.

Data assimilation (DA) offers a convenient framework to reduce model uncertainties by combining observations with the model simulation taking into account errors in both sources of information. Prior to DA, the main sources of uncertainties should be identified and included in the control vector; this is achieved with a sensitivity analysis (SA) study that allows to classify uncertainties in the inputs with respect to their impact on the model outputs, for instance in terms of variance using Sobol’ indices (Iooss and Saltelli 2016). Several studies in the framework of hydraulics demonstrated the merits of DA (Barthélémy 2015; Cloke and Pappenberger 2009; Dechant and Moradkhani 2011; Habert et al. 2016; Moradkhani et al. 2005) to provide a more accurate river state. This is achieved by inferring an optimal set of parameters (e.g. river and floodplain friction coefficients, upstream and lateral river discharge, bathymetry) and/or by simulating a more accurate river state, thus constituting a paradigm shift for real-time flood forecasting. Ensemble-based methods such as the Ensemble Kalman Filter (EnKF) (Durand et al. 2008; Moradkhani et al. 2005) and the Particle Filter (PF) (Matgen et al. 2010; Parrish et al. 2012) are popular algorithms; they articulate as a two-step procedure derived from Bayesian inference: (1) a forecast step to sample the uncertain inputs and propagate the uncertainty through the model, thus providing an ensemble of river states; and (2) an analysis step to weight each member or particle of the ensemble based on its discrepancies to the available observations and to derive in the case of parameter estimation, a correction on the inputs that is then propagated to the river states by model integration. In the EnKF algorithm, the weights are provided by the stochastic estimation of covariance matrices between errors in model inputs and outputs. In contrast, in the PF algorithm, the weights correspond to likelihoods associated with the probability density function (PDF) of the control vector conditioned upon the observations; PF provides an alternative to EnKF when the model is subject to strong non-linearity and non-Gaussian errors.

Monte Carlo (MC) techniques are generic, robust and easily portable on massively parallel supercomputers; yet they remain computationally expensive due to their slow convergence rate scaling as the inverse of the square root of the number of particles (Li and Xiu 2008). As shown in Barthélémy (2015) and Bozzi et al. (2014), a large number of forward model evaluations should be carried out to converge the stochastic evaluation of error statistics such as PDFs, Sobol’ indices and covariance matrices. There is therefore a need to develop efficient and robust uncertainty quantification (UQ) methods in the context of DA for hydraulics to limit (1) the number of significant sources of uncertainties and (2) the computational cost of quantifying uncertainties on the river state, e.g. moments (mean, covariance) and PDF, while preserving the accuracy of the mapping $\mathcal{M}$ between the uncertain inputs $\mathbf{x}$ and the vector of $M$ river water heights $\mathbf{h}$:

[TABLE]

The key idea of non-intrusive UQ methods is to build a cost-effective surrogate, also called reduced model, metamodel or response surface, mimicking the mapping $\mathcal{M}$ to perform UQ and SA steps (Iooss et al. 2010; Iooss and Saltelli 2016; Lamboni et al. 2011; Le Maitre and Knio 2010; Saltelli et al. 2007; Storlie et al. 2009). Formulating the surrogate model relies on a limited number of forward model integrations referred to as the training sample $(\mathcal{X},\mathcal{Y})=\left(\mathbf{x}^{(k)},\mathbf{h}^{(k)}\right)_{1\leq k\leq N}$ where $N$ is the training sample size and $\mathbf{h}^{(k)}:=\mathcal{M}(\mathbf{x}^{(k)})$ corresponds to the deterministic integration of the forward model $\mathcal{M}$ as a black box for the $k$th set of input parameters $\mathbf{x}^{(k)}$. Several surrogate models are found in the literature, among whom generalized linear models, polynomial models, splines, polynomial chaos expansions, artificial neural networks, Gaussian process models, …

Polynomial chaos (PC) approach has received much attention lately (Dubreuil et al. 2014; Sudret 2008; Xiu 2010; Xiu and Karniadakis 2002). The PC surrogate model is formulated as a polynomial expansion, in which the basis is defined according to the distribution of the uncertain inputs $\mathbf{x}$ and the coefficients relate to the statistics of the output $\mathbf{h}$. The coefficients are computed so as to fit the training set $(\mathcal{X},\mathcal{Y})$, either using regression or spectral projection methods. The merits of PC surrogates were demonstrated in various fields, e.g. structural mechanics (Dubreuil et al. 2014; Berveiller 2005), computational fluid dynamics (Hosder and Walters 2006; Lucor et al. 2007; Saad 2007), hydrology (Deman et al. 2015), hydraulics (Liang et al. 2008), wildfires (Rochoux et al. 2014). A complementary approach between PC and EnKF was presented in Li and Xiu (2009) and tested in the framework of wildfire behaviour forecasting in Rochoux et al. (2014). A PC surrogate was used in place of the forward model to estimate the Kalman gain and thereby design a cost-effective EnKF inferring new estimates of the input parameters (e.g. biomass moisture and fuel aspect ratio) reducing bias in the fire front location prediction. The resulting algorithm achieved convergence in spite of the non-linear surface response, showing promising results for environmental risk monitoring.

Gaussian processes (GP) that are strongly related to Kriging in geostatistics also are of increasing interest (Le Gratiet et al. 2017; Lockwood and Anitescu 2012; Rasmussen and Williams 2006). The GP formalism treats the forward model response as a realization of a Gaussian stochastic process indexed by $\mathbf{x}$ and fully characterized with mean and covariance functions conditioned by the training set $(\mathcal{X},\mathcal{Y})$. The GP surrogate is built first, by defining a covariance kernel function between output values and a trend function and then, by estimating the hyperparameters (e.g. variance, characteristic length scale) that provide a good fit to the training set. GP surrogates were introduced in the context of SA for estimating Sobol’ indices (Oakley and O’Hagan 2004; Marrel et al. 2009; Le Gratiet et al. 2014).

PC and GP surrogate models have recently been compared for UQ and SA studies (Le Gratiet et al. 2017; Owen et al. 2015; Schoebi et al. 2015). Owen et al. (2015) evaluated the performance of each type of surrogate in terms of output mean, variance and PDF estimation. Le Gratiet et al. (2017) compared Sobol’ indices with applications in structural mechanics; PC and GP surrogate models were found to feature a similar accuracy with respect to the $Q_{2}$ validation coefficient for a given size of the training set. They also emphasized that the ranking between PC and GP approaches remains problem-dependent. Since the application of PC surrogate models to hyperbolic conservation laws such as SWE remains a challenging task (Després et al. 2013; Birolleau et al. 2014), it is of great interest to evaluate the performance of PC and GP approaches when applied to open-channel flow simulation in fluvial or critical configurations.

In this paper, our objective is to evaluate the performance of PC and GP surrogates mimicking the behaviour of 1-D open-channel flows for UQ and SA steps that are important in the design of EnKF and PF algorithms. The quantity of interest is the 1-D water level discretized along 50 km of the Garonne River (South-West France). As a preliminary step, the focus is made on steady flows and the flow remains fluvial. The treatment of strong non-linearity induces by fluvial/critical transitions is beyond the scope of this study. The main sources of uncertainties are the upstream discharge that is constant in steady state conditions and the friction coefficient that is a piecewise constant function. Under these hypotheses, the present work features scalar inputs and spatially varying outputs due to the heterogeneous bathymetry of the Garonne River. A convergence study is first carried out to determine the size of the training set $N$ that is required to build a valid surrogate using either PC expansion or GP model. GP and PC surrogates are then compared when a computational budget is set, i.e. for the same number of snapshots (also called simulations in the paper) $N$ used to construct each surrogate. The comparison is carried out with respect to the following metrics on $\mathbf{h}$: PDF that is of interest for Bayesian inference and PF algorithms; spatially varying Sobol’ indices (associated with the correlation between each uncertain input and the spatially distributed output) and correlation matrix (associated with the spatial correlation of the output) that are of interest for variational and EnKF algorithms.

The structure of the paper is as follows. Section 2 introduces the 1-D SWE solved using MASCARET and the Garonne River case study. The PC and GP techniques are detailed in Sect. 3 along with error metrics. Section 4 presents the results of the comparative study between PC and GP with respect to a MC reference. Conclusions and perspectives are given in Sect. 5.

2 Hydraulic Modelling

The SWE represent the dynamics of open-channel flows, typically in rivers with small bathymetry variations (Horritt and Bates 2002). They form a hyperbolic system of partial differential equations that characterize subcritical and supercritical flows subject to hydraulic jumps. Here, we only deal with subcritical flows, in the plain in the downstream part of the Garonne River.

2.1 1-D shallow water equations (SWE)

We consider a 1-D hydraulic model commonly used in hydraulic engineering and flood forecasting. The stream channel is described by a hydraulic axis corresponding to the main flow direction, implying that the river channel is represented by a series of cross-sections (or profiles) identified by a curvilinear abscissa $a$ ranging from $a_{\text{in}}$ upstream of the river to $a_{out}$ downstream. 1-D SWE are derived from mass conservation and momentum conservation. The equations are written in terms of discharge (or flow rate) $Q$ (m3 s*-1*) and hydraulic section $A$ (m2) that relates to water level (or water height) $h$ (m) such that $A\equiv A(h)$. The non-conservative form of the 1-D SWE for non-stationary flow reads (Thual 2010)

[TABLE]

with $g$ the gravity, $S_{0}$ the channel slope and $S_{f}$ the friction slope. In the present study, the SWE are combined with the Manning-Strickler formula to parameterize the friction slope $S_{f}$ such as:

[TABLE]

where $R(h)=A(h)/P(h)$ (m) is the hydraulic radius written as a function of the wet perimeter $P(h)$, and where $K_{s}$ (m1/3 s${}^{-1})$ is the Strickler friction coefficient. The pair $(h,Q)$ forms the hydraulic state varying in time and space. For steady flows, Eq. 4 simplifies to:

[TABLE]

where $Fr$ is the dimensionless Froude number

[TABLE]

The smooth solutions for Eq. 9 are called backwater curves when the downstream boundary condition is prescribed in a deterministic way. To solve Eq. (9), the following input data are required: bathymetry, upstream and downstream boundary conditions, lateral inflows, roughness coefficients and initial condition for the hydraulic state. The imperfect description of this data translates into errors in the simulated hydraulic state $(h,Q)$. To understand the structure of these errors, it is of prime importance to determine which input variables contribute, and to what extent, to the variability in the hydraulic state at different curvilinear abscissas $a$ along the river channel, for instance via a SA study (Sect. 3.3).

We use the MASCARET software to simulate the 1-D SWE in Eq. (4) and predict $(h,Q)$ along the discretized curvilinear abscissa of the hydraulic network $a\in[a_{\text{in}},a_{\text{out}}]$. The SWE are solved here with the steady kernel of MASCARET based on a finite difference scheme (Goutal and Maurel 2002; Goutal et al. 2012), meaning that $(h,Q)$ only varies in space. MASCARET is part of the TELEMAC-MASCARET open-source modelling package developed at EDF (Électricité de France R&D) in collaboration with CEREMA (Centre d’Étude et d’expertise sur les Risques, l’Environnement, la Mobilité et l’Aménagement); it is commonly used for dam-break wave simulation, reservoirs flushing and flooding.

2.2 Garonne River Case Study

The present study is carried out on a real hydraulic network over the Garonne River in South-West France. The Garonne river flows along 647 km from the Pyrenees to the Atlantic Ocean draining an area of $55,000$ km2 (corresponding to the fourth-largest river in France). The present study focuses on a 50 km reach from Tonneins ($a_{\text{in}}=13$ km) to La Réole ($a_{\text{out}}=62$ km) with one observing station at Marmande ($a=36$ km) as presented in Fig. 1a. The mean slope over the reach is $S_{0}=3.3$ m km*-1*; the mean width of the river is $W=250$ m; the bank-full discharge is approximately $\overline{Q}=2550$ m3 s*-1*. Despite the existence of active floodplains, this reach can be accurately modelled by a 1-D hydraulic model using Eq. (4) (Besnard and Goutal 2011).

Figure 1b presents the non-uniform bathymetry profile along the 50 km reach, interpolated from $83$ on-site bathymetry cross-sections. Friction for the river channel and its floodplain is prescribed over three zones separated by dashed lines. The Strickler coefficients $K_{s_{1}}$, $K_{s_{2}}$ and $K_{s_{3}}$ are used to characterize friction through Eq. (5) and are uniform per zone. The observing station at Marmande is located at the beginning of the third zone associated with $K_{s_{3}}$.

The upstream steady boundary condition is prescribed by $Q(a_{\text{in}})=Q_{\text{in}}$; the discharge $Q$ is constant along the reach ($Q=Q_{\text{in}}$). The downstream boundary condition is prescribed with a local rating curve $RC$ established at La Réole that sets $h(a_{\text{out}})=RC(Q_{\text{out}})=h_{\text{out}}$. The hydraulic model has been calibrated using channel and floodplain roughness coefficients as free parameters (Besnard and Goutal 2011).

2.3 Sources of uncertainties and Quantity of Interest

UQ and DA for flood forecasting are essential to ensure the predictive capability of the surrogate model at the observing stations such as Marmande. Previous work (El Moçayd et al. ) has shown that in steady state conditions and given some assumptions on the statistics of the input uncertain variables, the sensitivity of the hydraulic state at Marmande to $K_{s_{3}}$ is predominant; the sensitivity at Marmande to $K_{s_{1}}$ and $K_{s_{2}}$ is negligible. Hence, the main sources of uncertainties taken into account here are the upstream mass flow rate $Q$ and the Strickler coefficient $K_{s_{3}}$. We denote by $\mathbf{x}=(Q,K_{s_{3}})$ the random vector of size $d=2$. From expert knowledge, $Q$ and $K_{s_{3}}$ are considered as independent random variables where $Q$ (m3 s*-1*) follows the normal distribution $\mathcal{N}(4031,400)$ and $K_{s_{3}}$ (m1/3 s*-1*) follows the uniform distribution $\mathcal{U}(15,60)$.

MASCARET provides as output the water height over 463 cross-sections for the Garonne case. In this work, we focus on the water height at $M=14$ stations evenly distributed along the 50 km reach, among which Marmande at $a=36$ km. We denote by $\mathbf{h}$ the vector of $M$ observed water levels for one realization of MASCARET. Note that the UQ methods we present in Sec. 3 are highlighted in this work for $M=14$ but they can easily be applied to a case where the size of the state vector $M$ is much larger than the size of the training set $N$.

A database noted $\mathcal{D}_{N_{\text{ref}}}$ and containing $N_{\text{ref}}=\numprint{100,000}$ MASCARET simulations was compiled as a reference for the study. Each simulation corresponds to a different pair of inputs ($Q,K_{s_{3}}$) resulting from a MC random sampling following $Q\sim\mathcal{N}(4031,400)$ and $K_{s_{3}}\sim\mathcal{U}(15,60)$. In the present study $D_{N_{\text{ref}}}$ is partly used to build the PC and GP surrogates via the training set $\mathcal{D}_{N}=(\mathcal{X}$, $\mathcal{Y})\subset\mathcal{D}_{N_{\text{ref}}}$ of size $N$; it is fully used to validate them a posteriori over a large ensemble (Sec. 3.4).

3 Surrogate Models

Two types of surrogate models are used in the present work to approximate the behaviour of MASCARET with respect to $\mathbf{x}=(Q,K_{s_{3}})$ (the size of the input uncertainty space is $d=2$): PC expansion on the one hand, GP regression combined with Proper Orthogonal Decomposition (POD) further denoted by pGP on the other hand. The common idea of PC and pGP approaches is to design for each curvilinear abscissa $a\in\{a_{1},\cdots,a_{M}\}$ a surrogate water level with a weighted finite sum of basis functions:

[TABLE]

where the coefficients $\gamma_{a,i}$ and the basis functions $\Psi_{i}$ are calibrated by the training set $\mathcal{D}_{N}$. For PC expansion, the maximum order of the decomposition and the basis functions are first chosen, then a spectral projection strategy is implemented to compute the decomposition coefficients in the polynomial basis. The pGP strategy is two-fold: the first step consists in transforming the sampled output space to an orthogonal one, eventually reducing its dimension; the second step consists in interpolating with a GP approach the principal components of this new space to express water level in the basis for any friction and discharge values. It should be noted that the GP model could be replaced by any other interpolator. The main difference between PC (Sect. 3.1) and pGP (Sect. 3.2) approaches stands in the nature of these models: pGP interpolates the training points and captures local variations, while PC is a regression method focusing on the global behaviour of the model. Basis functions and calibration methods also differ between these two approaches.

3.1 Polynomial Chaos (PC) surrogate model

The algorithm to build the PC surrogate proceeds as follows:

choose the polynomial basis $\{\Psi_{i}\}_{i\geq 0}$ according to the assumed PDF of the inputs $\mathbf{x}=(Q,K_{s_{3}})$, 2. 2.

choose the total polynomial degree $P$ according to the complexity of the physical processes, 3. 3.

truncate the expansion to $r_{\text{pc}}$ terms to keep the predominant information given by the forward model using standard truncation strategy ($r_{\text{pc}}$ depends on $d$ and $P$), 4. 4.

apply spectral projection strategy (i.e. Gaussian quadrature rule) to compute the coefficients $\{\gamma_{a,i}\}_{i\in\mathbb{N}^{d}\atop|i|\leq P}$ for each curvilinear abscissa $a$ using $N=(P+1)^{d}$ snapshots from the simulation database $\mathcal{D}_{N_{\text{ref}}}$, 5. 5.

formulate the surrogate model $\mathcal{M}_{\text{pc}}$ at each curvilinear abscissa $a$, which can be evaluated for any new pair of parameters $\mathbf{x}^{*}=(Q^{*},K_{s_{3}}^{*})$.

Note that we use standard truncation and projection strategies presented in Le Maitre and Knio (2010) and Xiu (2010).

3.1.1 Polynomial basis

Each component of the random vector $\mathbf{x}$ defined in the input physical space is standardized and noted $\boldsymbol{\zeta}$ in the following way: $\zeta_{i}=\frac{x_{i}-\mu_{i}}{\sigma_{i}}$ where $\mu_{i}=N^{-1}\sum_{k=1}^{N}x_{i}^{(k)}$ and $\sigma_{i}=\sqrt{(N-1)^{-1}\sum_{k=1}^{N}\left(x_{i}^{(k)}-\mu_{i}\right)^{2}}$. Assuming that the SWE solution is of finite variance, each component $h_{a}$ of the water level vector $\mathbf{h}$ can be considered as a random variable for which there exists a polynomial expansion of the form (11) that represents how the water level $h_{a}$ varies according to changes in $Q$ and $K_{s_{3}}$.

$h_{a}$ is projected onto a stochastic space spanned by the orthonormal polynomial functions $\{\Psi_{i}\}_{i\geq 0}$. These functions are orthonormal with respect to the joint density $\rho(\boldsymbol{\zeta})$, i.e.

[TABLE]

with $\delta_{ij}$ the Kronecker delta function and $Z\subseteq\mathbb{R}^{d}$ the space in which $\boldsymbol{\zeta}$ evolves. In practice, the orthonormal basis is built using the tensor product of 1-D polynomial functions: $\Psi_{i}=\Psi_{i,1}\ldots\Psi_{i,d}$ where $i$ is the multi-index $(i_{1},\ldots,i_{d})\in\{0,1,\cdots,P\}^{d}$. The choice for the basis functions depends on the probability measure of the random variables. According to Askey’s scheme, the Hermite polynomials form the optimal basis for random variables following the standard Gaussian distribution, and the Legendre polynomials are the counterpart for the standard uniform distribution (Xiu and Karniadakis 2002).

3.1.2 Truncation Strategy

In practice, the sum in Eq. (11) is truncated to a finite number of terms $r_{\text{pc}}$. Using a standard truncation strategy $r_{\text{pc}}$ is constrained by the number of random variables $d$ and by the total polynomial degree $P$ as:

[TABLE]

meaning that all polynomials involving the $d$ random variables of total degree less or equal to $P$ are retained in the PC expansion. The PC approximated water level at curvilinear abscissa $h_{\text{pc}}(a)$ is formulated as:

[TABLE]

Note that for small $d$, advanced truncation strategies that consist in eliminating high-order interaction terms or using sparse structure (Blatman 2009; Migliorati et al. 2013) are not necessary.

3.1.3 Spectral projection strategy

We focus here on non-intrusive approaches to numerically compute the coefficients $\{\gamma_{a,i}\}_{i\in\mathbb{N}^{d}\atop|i|<P}$ in Eq. (14) using $N$ snapshots from $\mathcal{D}_{N}$. The spectral projection relies on the orthonormality property of the polynomial basis. The $i$th coefficient $\gamma_{a,i}$ is computed using Gaussian quadrature as:

[TABLE]

where $\mathbf{h}^{(k)}=\mathcal{M}(\mathbf{x}^{(k)})$ is the snapshot $\mathcal{D}_{N}$ corresponding to the MASCARET simulation for the $k$th quadrature root $\mathbf{x}^{(k)}$ of $\Psi_{i}$ (in the physical space), and where $w^{k}$ is the weight associated with $\mathbf{x}^{(k)}$. $(P+1)$ is the number of quadrature roots required in each uncertain direction to ensure an accurate calculation of the integral $<h_{a},\Psi_{i}>$. Hence, $N=(P+1)^{2}$ for the PC surrogates built in this study.

3.2 POD-based Gaussian Process (pGP) surrogate

The algorithm to build a Gaussian Process surrogate relies on a POD and proceeds as follows:

choose the size of the training set $N$, 2. 2.

draw $N$ samples (or snapshots) in the input random space $\mathbf{x}=(Q,K_{s_{3}})$ with Halton’s low discrepancy sequence111From a Halton sequence, the nearest points (standardized Euclidean distance) of the database $\mathcal{D}_{N}$ are chosen. from the simulation database $\mathcal{D}_{N}$, 3. 3.

formulate the centred snapshot matrix $\mathbf{Y}$ from the $N$ water level snapshots, 4. 4.

achieve a POD on $\mathbf{Y}$ using the snapshot method to derive the basis vectors $\{\Psi_{i}\}$ and the corresponding coefficients $\{\gamma_{a,i}\}$ (any snapshot can be expressed as a linear combination of the basis vectors and coefficients), 5. 5.

replace each basis vector $\{\Psi_{i}\}$ associated to the $N$ snapshots by a basis function via Gaussian Process regression for any $\mathbf{x}^{*}$, 6. 6.

formulate the surrogate model $\mathcal{M}_{\text{gp}}$ for water level at each curvilinear abscissa $a$, which can be evaluated for any $\mathbf{x}^{*}$.

Note that we follow the choices made by Braconnier et al. (2011).

3.2.1 Snapshot method

The key idea of the snapshot method (Sirovich 1987) is to achieve a POD of the centred snapshot matrix $\mathbf{Y}\in\mathbb{M}_{M,N}(\mathbb{R})$, which gathers the water level computed at each curvilinear abscissa for the $N$ snapshots, from which the sample mean is subtracted. For simplicity purpose, the water level anomaly at curvilinear abscissa $a$ is denoted $h$ in the following. Thus, $\mathbf{Y}=\left(h_{a_{i}}^{(j)}\right)_{1\leq i\leq M\atop 1\leq j\leq N}$. The snapshots correspond to the column vectors; the $k$th snapshot of size $M=14$ is denoted by $\mathbf{h}^{(k)}$.

Based on many observations of a random vector, the POD gives the orthogonal directions of largest variances (or modes) in the probabilistic vector space in order to reduce the vector space dimension (Chatterjee, Anindya 2000). Note that for simplicity purpose, the adjective centred is dropped in the following when referring to the centred snapshot matrix $\mathbf{Y}$.

The POD of the snapshot covariance matrix $\mathbf{C}=N^{-1}\,\mathbf{Y}^{\mathrm{T}}\,\mathbf{Y}\in\mathbb{M}_{N}(\mathbb{R})$ is equivalent to the Singular Value Decomposition (SVD) of the snapshot matrix $\mathbf{Y}$:

[TABLE]

where $\mathbf{U}\in\mathbb{M}_{M}(\mathbb{R})$ is an orthogonal matrix diagonalizing $\mathbf{Y}\mathbf{Y}^{\mathrm{T}}$ ($\mathbf{u}_{k}$, the $k$th column of $\mathbf{U}$, is a left singular vector of $\mathbf{Y}$), where $\mathbf{V}\in\mathbb{M}_{N}(\mathbb{R})$ is an orthogonal matrix diagonalizing $\mathbf{Y}^{\mathrm{T}}\mathbf{Y}$ ($\mathbf{v}_{k}$, the $k$th column of $\mathbf{V}$, is a right singular vector of $\mathbf{Y}$), and where $\mathbf{\Lambda}\in\mathbb{M}_{M,N}(\mathbb{R})$ is a rectangular diagonal matrix including $r_{p}=\min(M,N)$ singular values on its diagonal. The singular values $\{\lambda_{k}\}_{1\leq k\leq r_{p}}$ are the square roots of the eigenvalues of $\mathbf{C}$. Note that in this study, we do not reduce further the rank of the snapshot matrix $\mathbf{Y}$. Since the number of stations ($M=14$) is lower than the size of the training set $N$, the rank of $\mathbf{Y}$ is here $r_{p}=M$.

For a given curvilinear abscissa $a$, any snapshot $h_{a}(\mathbf{x}^{(k)})$ can then be retrieved as a linear combination of $r_{p}$ modes $\{\Psi_{i}\}_{1\leq i\leq r_{p}}$:

[TABLE]

where for any $i\in\{1,\ldots,M\}$, $\gamma_{a,i}:=U_{a,i}$ and $\Psi_{i}(\mathbf{x}^{(k)}):=(\boldsymbol{\Lambda}\mathbf{V}^{T})_{i,k}$.

In the present study, the size of the state vector $M$ is smaller than the size of the training set $N$, implying that there are $r_{p}=M$ non-zero singular values. This is usually not the case, but this does not change the methodology. As the objective is to compare PC and pGP and as there is a limited number of stations $M$, all modes are kept to avoid the loss of information, there is no need to reduce the rank of the singular matrix. Note that when moving to 2-D or 3-D cases (beyond the scope of the present work), the size of the state vector could be of the order of thousands components, making dimension reduction necessary.

3.2.2 Regression Procedure

Our objective is now to generalize the data set $\{\Psi_{i}(\mathbf{x}^{(k)})\}_{1\leq k\leq N}$ associated with the design of experiments $\mathcal{X}$ to any new input vector $\mathbf{x}^{*}$, in order to approximate the model output $h_{a}$ at any curvilinear abscissa $a$. We thus propose the following surrogate model based on GP regression and Eq. (17):

[TABLE]

where $\Psi_{\text{gp},i}$ is a GP model calibrated from the training set $\left\{\mathbf{x}^{(k)},\Psi_{i}(\mathbf{x}^{(k)})\right\}_{1\leq k\leq N}$ as detailed in the following.

As stated by Rasmussen and Williams (2006), a GP is a random process (here the mode $\Psi_{i}$) indexed over a domain (here $\mathbb{R}^{d}$), for which any finite collection of process values (here $\left\{\Psi_{i}(\mathbf{x}^{(k)})\right\}_{1\leq k\leq N}$) has a joint Gaussian distribution. Concretely, let $\widetilde{\Psi}_{i}$ be a Gaussian random process fully described by its zero mean and its correlation $\pi_{i}$:

[TABLE]

with $\pi_{i}(\mathbf{x},\mathbf{x}^{\prime})=\mathbb{E}\left[\widetilde{\Psi}_{i}(\mathbf{x})\widetilde{\Psi}_{i}(\mathbf{x}^{\prime})\right]$. In our case, the correlation function $\pi$ (or kernel) is chosen as a squared exponential:

[TABLE]

where $\ell_{i}$ is a length scale describing dependencies of model output between two input vectors $\mathbf{x}$ and $\mathbf{x}^{\prime}$, and where $\sigma_{i}^{2}$ is the variance of the output signal. Then the surrogate model of interest is the mean of the GP resulting of conditioning $\widetilde{\Psi}_{i}$ by the training set $\left\{\Psi_{i}\left(\mathbf{x}^{(k)}\right)\right\}_{1\leq k\leq N}$. For any $\mathbf{x}^{*}\in\mathbb{R}^{d}$,

[TABLE]

where $\beta_{k,i}=\left(\mathbf{\Pi}_{i}+\tau^{2}\,\mathbf{I}_{N}\right)^{-1}\left(\Psi_{i}\left(\mathbf{x}^{(1)}\right)\ldots\Psi_{i}\left(\mathbf{x}^{(N)}\right)\right)^{T}$ with $\mathbf{\Pi}_{i}=\left(\pi_{i}\left(\mathbf{x}^{(j)},\mathbf{x}^{(k)}\right)\right)_{1\leq j,k\leq N}$, and where $\tau$ (referred to as the nugget effect) avoids ill-conditioning issues for the matrix $\mathbf{\Pi}$. The hyperparameters $\left\{\ell_{i},\sigma_{i},\tau\right\}$ are optimized by maximum likelihood applied to the data set $\mathcal{D}_{N}$ using a basin hopping technique (Wales and Doye 1997).

3.3 Statistical Analysis

In the following, the water level at curvilinear abscissa $a$ estimated with either the PC surrogate or the pGP surrogate is denoted for simplicity purpose by $\widehat{h}_{a}$, while $h_{a}$ denotes the water level obtained using the MASCARET forward model.

3.3.1 Statistical moments and PDF

Using a standard MC approach on the reference sampling $\mathcal{D}_{N_{\text{ref}}}$, the water level mean value and standard deviation (STD) at a given curvilinear abscissa $a$ noted $(\mu_{h_{a}},\sigma_{h_{a}})$ are formulated as:

[TABLE]

The covariance matrix of the simulated water level denoted by $\mathbf{C}\in\mathcal{M}_{M}(\mathbb{R})$ is stochastically formulated as:

[TABLE]

with $\mathbf{h}^{(k)}=\mathcal{M}(\mathbf{x}^{(k)})$ the $k$th sample in $\mathcal{D}_{N_{\text{ref}}}$ containing the water level at the $M$ observed curvilinear abscissas and with $\overline{\mathbf{h}}=N_{\text{ref}}^{-1}\sum_{k=1}^{N_{\text{ref}}}\mathbf{h}^{(k)}$ the ensemble mean.

Using the PC surrogate method, the statistical moments can be derived analytically from the coefficients $\{\gamma_{a,i}\}_{i\in\mathbb{N}^{d}\atop|i|<P}$ such that at a given curvilinear abscissa $a$, the water level mean and STD ($\mu_{h_{pc,a}}$, $\sigma_{h_{\text{pc},a}}$) read:

[TABLE]

The covariance matrix $\mathbf{C}_{\text{pc}}$ can also be directly computed from the coefficients:

[TABLE]

with $\text{cov}\left(h_{pc,a_{m}},h_{pc,a_{n}}\right)$ the matrix component corresponding to the water level covariance at grid points $a_{m}$ and $a_{n}$.

Using the pGP surrogate method, both statistical moments and covariances are stochastically estimated as in the MC approach, meaning that $h$ is replaced by the surrogate $h_{\text{pgp}}$ in Eqs. (22)–(24).

For the MC random sampling method as well as for the PC and pGP surrogate methods, the PDF of the water level at each of the $M=14$ stations is reconstructed using a kernel smoothing procedure on a large enough stochastic sampling of $\mathcal{M}$, $\mathcal{M}_{\text{pc}}$ or $\mathcal{M}_{\text{gp}}$, respectively (Wand and Jones 1995; Hastie et al. 2009).

3.3.2 Sensitivity Analysis

Sobol’ indices (Sobol 1993; Saltelli et al. 2007) are commonly used for SA based on variance analysis. They provide the quantification of how much of the variance of the quantity of interest is due to the uncertainty in the input parameters assuming these random variables are independent and the random output is squared integrable. The water level (simulated either by MASCARET $h_{a}$ or by the PC or pGP surrogate $\widehat{h}_{a}$) variance at curvilinear abscissa $a$ decomposes as:

[TABLE]

where $\mathbb{V}_{i}(\widehat{h}_{a})=\mathbb{V}\left[\mathbb{E}(\widehat{h}_{a}|x_{i})\right]$, $\mathbb{V}_{ij}(\widehat{h}_{a})=\mathbb{V}\left[\mathbb{E}(\widehat{h}_{a}|x_{i}x_{j})\right]-\mathbb{V}_{i}(\widehat{h}_{a})-\mathbb{V}_{j}(\widehat{h}_{a})$ and more generally, for any $I\subset\{1,\ldots,d\}$, $\mathbb{V}_{I}(\widehat{h}_{a})=\mathbb{V}\left[\mathbb{E}(\widehat{h}_{a}|x_{I})\right]-\sum_{J\subset I\,\text{s.t.}\,J\neq I}\mathbb{V}_{J}(\widehat{h}_{a})$. The Sobol’ indices for water level at curvilinear abscissa $a$ with respect to the $i$th and $j$th components of the random input vector read:

[TABLE]

$S_{i}^{a}$ is the first order Sobol’ index corresponding to the ratio of output variance due to the $i$th input parameter uniquely, and $S_{ij}$ is the second-order Sobol’ index describing the ratio of output variance due to the $i$th parameter in interaction with the $j$th parameter. Also the total Sobol’ index that corresponds to the whole contribution of the $i$th input parameter reads:

[TABLE]

For both MC and pGP approaches, the Sobol’ indices are stochastically estimated. The conditional variances are estimated using Martinez’ formulation. This estimator is stable and it provides asymptotic confidence intervals for first order and total order indices (Baudin et al. 2016).

For the PC approach, Sobol’ indices can be directly derived from the PC coefficients. For the $i$th component of the input random variable $\mathbf{x}$, the Sobol’ index reads:

[TABLE]

with $\sigma_{h_{\text{pc},a}}$ the STD computed in Eq. (26).

3.4 Error Metrics

In the present study, three error metrics are used to assess the quality of the surrogate water level at a given curvilinear abscissa $\widehat{h}_{a}$: the $Q_{2}$ predictive coefficient (Marrel et al. 2009), the Kolmogorov-Smirnov test to evaluate the similarity between PDF (Clarke et al. 1992) and the Root Mean Square Error (RMSE). The validation is carried out over a full database $D_{N_{\text{ref}}}$ of size $N_{\text{ref}}$.

Predictive coefficient $Q_{2}$

At a given curvilinear abscissa $a$, the $Q_{2}$ predictive coefficient reads:

[TABLE]

Kolmogorov-Smirnov statistical test

Let $T_{F}$ (resp. $T_{G}$) be a random variable with cumulative distribution function (CDF) $F$ (resp. $G$). Let $F_{n}$ (resp. $G_{m}$) be its empirical CDF built from $n$ (resp. $m$) independent realizations of $T_{F}$ (resp. $T_{G}$). Then, let us define the test statistics:

[TABLE]

The null hypothesis for the Kolmogorov-Smirnov statistical test supposes that $T_{F}$ and $T_{G}$ are identically distributed, i.e. $F=G$. The Kolmogorov-Smirnov test leads us to reject this hypothesis with a type I error $\alpha\in]0,1[$ when:

[TABLE]

with $c(\alpha)$ a tabulated value found in the literature (Smirnov 1939). In the present study, using $\alpha=0.05$ and $n=m=N_{\text{ref}}$, the null hypothesis is rejected if $D>\numprint{6.082e-3}$.

Root mean square error (RMSE)

The RMSE is used to evaluate the accuracy in the estimation of the correlation matrix and of the Sobol’ indices with respect to the MC stochastic estimation over the data set $D_{N_{\text{ref}}}$, i.e.

[TABLE]

where $\Box$ designates one component of the correlation matrix or one Sobol’ index at the given curvilinear abscissa $a$.

3.5 Numerical implementation and validation

Due to the increasing interest in UQ over the last decade, a significant number of UQ-dedicated tools/libraries are now available for the scientific community. OpenTURNS (see www.openturns.org) is an open-source (GNU LGPL) scientific library developed since 2005 by EDF, Airbus, Phimeca and IMACS and usable as a Python module dedicated to uncertainty treatment and reliability analysis in a structured industrial approach (Baudin et al. 2015). OpenTURNS offers a wide catalogue of features for UQ (e.g. PC expansion, Dutka-Malen et al. 2009) and benefits from a well-organized developers’ and users’ community (forum, training, user guides). It can either be used as a Python module or as a component within a coupling platform. For further information on the OpenTURNS library, the reader is referred to the online reference and use case guides.

In the present study, MASCARET and OpenTURNS are both integrated components of the SALOME platform developed at EDF (see www.salome-platform.org). This integrated framework allows for an efficient use of MASCARET as a Python function for simulating an ensemble of water levels and for building the PC surrogate models using the methods implemented in OpenTURNS. The pGP surrogate models rely on the use of the JPOD Python tool developed at CERFACS (Braconnier et al. 2011) since 2007 and recently augmented with OpenTURNS’ UQ capabilities (Roy 2016). It acts as a Python platform that is able to interact with an external simulation code, i.e. generate a design of experiment, run the corresponding simulations and perform statistical analysis of the outputs. The GP implementation relies on the package Scikit-learn (Pedregosa et al. 2012).

Both PC and pGP strategies were validated on: (i) classical optimization functions with scalar output such as Ishigami (Ishigami and Homma 1990) and Michalewicz (Molga and Smutnicki 2005); and (ii) for a functional output case using simplified 1-D open-channel flow equations (backwater curves) (El Moçayd et al. ; Roy 2016). JPOD has been successfully tested and validated on industrial configurations during the European SimSAC and ALEF (FP7) projects.

4 Comparison of Polynomial Chaos (PC) and Gaussian Process (pGP) surrogates

4.1 Reference Monte Carlo (MC) results

We first present the results obtained using the MC reference sampling ($N_{\text{ref}}=100,000$ in $\mathcal{D}_{N_{\text{ref}}}$) in terms of water level mean, STD, PDF, correlation matrix as well as Sobol’ indices associated with $Q$ and $K_{s_{3}}$. These results are used as reference to evaluate the accuracy of the PC and pGP surrogate models.

Figure 2a displays the water level PDF computed from $\mathcal{D}_{N_{\text{ref}}}$ data set integrated with a MC approach for the $M=14$ stations along the curvilinear abscissa $a\in[a_{\text{in}},a_{\text{out}}]$. The mean water level (Eq. 22) is represented with a thick black line; the interval between two STD (Eq. 23) is represented with dotted lines; and the minimum and maximum water level values are represented with dashed lines. The upstream part of the river is under the influence of the upstream forcing, the spread of the ensemble tends to increase. The downstream part of the river is under the influence of the downstream boundary condition where the water level and discharge are related by the local rating curve $RC$, the spread of the ensemble tends to decrease near the downstream boundary condition. At Marmande ($a=36$ km), the flow is complex due to strong variation of the local bathymetry (Fig. 1b), the ensemble spread is larger and the PDF plotted in Fig. 2b features two main modes due to the change in backwater curves solutions for subcritical flow.

The Sobol’ indices for the water level $S_{Q}$ and $S_{K_{s_{3}}}$ computed from $\mathcal{D}_{N_{\text{ref}}}$ with respect to $Q$ and $K_{s_{3}}$ are presented in Fig. 3a along the curvilinear abscissa $a$. These indices confirm the previously mentioned spatial sensitivity. The water level variance is mostly explained by the upstream discharge variability for $a=[0;30~{}\text{km}]$. It is then mostly explained by the Strickler coefficient variance for $a=[30;60~{}\text{km}]$. Near the downstream boundary condition, the water level is related to the discharge by the local rating curve $RC$, the very last part of the network is thus under the influence of the discharge. First and total order indices are equal, meaning that there is no correlation between the errors in the input parameters.

Figure 3b displays the water level correlation matrix along the 50 km reach (derived using Eq. 24) that is estimated from $\mathcal{D}_{N_{\text{ref}}}$. The $n$th column of the matrix describes the water level error correlations between one given location on the channel $a_{n}$ and the rest of the channel $a_{m}$ with $a_{m}\in[a_{\text{in}},a_{\text{out}}]$. By definition, the correlation is equal to 1 on the diagonal, it decreases when the distance between $a_{n}$ and $a_{m}$ increases. We first analyse the correlation function for a point located upstream of the river ($a_{i}=15~{}\text{km}$) where $S_{Q}=0.9$ and $S_{K_{s_{3}}}=0.1$. Water level errors are strongly correlated in the upstream part of the river, which is under the influence of the upstream discharge boundary condition, where $S_{Q}$ is large. Errors between $a=15$ km and the rest of the river tends to decorrelate when the influence of $K_{s_{3}}$ increases (i.e. where $S_{K_{s_{3}}}$ is larger). We then analyse the correlation function for Marmande ($a=36~{}\text{km}$) where $S_{Q}=0.15$ and $S_{K_{s_{3}}}=0.85$. The correlation between water level errors at Marmande and the rest of the river is large in the vicinity of Marmande, where the influence of $K_{s_{3}}$ prevails. It then decreases for upstream and downstream locations that are under the influence of the upstream discharge and the downstream rating curve $RC$ (where $S_{Q}$ is large), respectively.

4.2 Convergence Analysis for Surrogates

We now use the same metrics as for the MC random sampling approach in Sec. 4.1 to evaluate the accuracy of both PC and pGP surrogates. The surrogate models are built using a training set $(\mathcal{X},\mathcal{Y})$ of size $N$ that is much smaller than that of the reference sample $N_{\text{ref}}$; they are then validated with respect to the reference MC results. The PC surrogate is built using a Gaussian quadrature rule with $N=(P+1)^{2}$ particles in the training set ($P$ is the total polynomial degree to be determined). The pGP approach is built using an approximate low-discrepancy Halton’s sequence of the same budget as for the PC approach. This is approximate in the sense that we consider the closest values to the standard Halton’s sequence that are part of the data set $\mathcal{D}_{N_{\text{ref}}}$. The sensitivity to the value of $P$ and thus to the size of the training set $N$ is investigated.

Both surrogates are computed with a fixed budget $N$ of 49 and 121 MASCARET evaluations. For PC, this value of $N$ corresponds respectively to $P=6$ and $P=10$. The pGP and PC response surfaces at Marmande are presented in Fig. 3. The design of experiment is represented by black dots. The colour map is evaluated by sampling each surrogate over the full data set $\mathcal{D}_{N_{\text{ref}}}$ made of $N_{\text{ref}}=100,000$ particles. It is found that the water level increases with increasing discharge $Q$ and decreasing Strickler coefficient $K_{s_{3}}$, consistently with MASCARET behaviour. Due to the quadrature rule, increasing the number of snapshots (from $P=6$ in Fig. 3c to $P=10$ in Fig. 3d) allows to build a higher order PC surrogate valid on a wider input range for $Q$ that is described by a Gaussian PDF. For $P=10$, some of the quadrature roots are outside of the MC sample and require additional MASCARET evaluations to build the PC surrogate $\mathcal{M}_{p}c$. Looking at the pGP design of experiments, the input space interval for $Q$ has been arbitrarily fixed to optimally represent the PDF. Since we consider a Gaussian distribution, its range has been bounded to $[3000;5000~{}$m^{3},s^{-1}\text{,}\mathrm{]}$$. Following Chebyshev’s theorem, this leads to a 90 % confidence interval.

The distribution of the particles in the design of experiments used by PC and pGP differs. This is done on purpose based on previous performance tests carried out with respect to the training set size $N$ and evaluating the surrogate accuracy, see (El Moçayd et al. ) and (Roy 2016), respectively. On the one hand, the design of experiments for the PC surrogate is constrained by the PDF of the uncertain inputs. We use here a quadrature-based PC since it was found to be more cost-effective than regression-based PC for a small dimensional problem ($d=2$) on the Garonne case (El Moçayd et al. ). On the other hand, using the approximate Halton’s sequence is known to be accurate for pGP surrogate (Damblin et al. 2013). This is indeed useful to cover the uncertain space without any bias and to have low discrepancy, meaning that most of the quantity of interest variance is captured and that a good $Q_{2}$ criterion is achieved. The choice of the design of experiment will be driven in future work according to the study objectives. For instance, DA usually relies on the assumption of Gaussian PDF for the input parameters, so the design of experiments can use this prior information. For risk analysis, threshold values are paramount and a design of experiment accounting for parameter space extrema is required.

The $Q_{2}$ error between the surrogate water level and the forward model water level for $D_{N_{\text{ref}}}$ averaged over the river is given in Table 1 for different sizes of training set $N$ varying between 49 and 256. The error remains below $\numprint{e-3}$, even for $N=49$ snapshots; it is only slightly improved when the number of snapshots $N$ increases to 256 and beyond.

The water level PDF estimated with the PC and the pGP surrogate models based on $49,121,256$ snapshots are compared in Fig. 5 at the curvilinear abscissa $a=15~{}\text{km}$ (near upstream boundary condition) and $a=36~{}\text{km}$ (Marmande). In the upstream part of the river, the PDF is uni-modal and is well represented with a small number of snapshots for both surrogates. On the contrary, at Marmande, the dynamics of the flow is more complex and the PDF is bimodal. Both PC and pGP surrogates are able to retrieve the overall shape of the PDF at $a=15~{}\text{km}$ and $a=36~{}\text{km}$. The shape is more accurate when the number of snapshots $N$ increases. This is quantified using a Kolmogorov-Smirnov statistical test, which measures the fit between the water level CDF computed from each surrogate model and that computed from the reference MC. Table 2 indicates that the null hypothesis (Eq. 34) is rejected for both surrogates computed with 49 snapshots (when $D>\numprint{6.082e-3}$) and accepted when at least 121 snapshots are used. For $N=49$, Fig. 5b–d show that the location of the first mode is shifted compared to MC reference. Increasing $N$ to 121 enables the PC approach to correctly locate this mode, while it enables the pGP approach to represent the second mode with an accurate amplitude, leading to an accepted null hypothesis. Each approach presents a particular limitation: the first mode is not well positioned by the pGP approach; the second mode is not captured by the PC approach. As for the tail of the PDF, it is correctly represented by the PC surrogate, while the pGP surrogate slightly oscillates around the shape of the reference MC PDF.

The RMSE for the Sobol’ indices for both pGP and PC with a fixed computational budget of $N=121$ simulations is of the order of $10^{-2}$ when computed over the 50 km reach. Figure 6 displays the squared error for the Sobol’ indices along the curvilinear abscissa (Eq. 35); the spatial pattern is similar for both surrogates and the squared error for each index is larger where the Sobol’ indices are larger. The RMSE for the correlation matrix for both pGP and PC with $N=121$ simulations are equal to $\text{RMSE}_{\text{pc}}=\numprint{3.67e-4}$ and $\text{RMSE}_{\text{gp}}=\numprint{4.59e-3}$. The spatial distribution of the squared error is plotted in Fig. 7. These results confirm the good behaviour of both PC and pGP surrogates with respect to MASCARET. For both Sobol’ indices and correlation matrices, the PC surrogate slightly outperforms pGP.

5 Discussion and Conclusions

The purpose of this work was to compare two popular strategies for building surrogate models, Polynomial Chaos (PC) and POD-based Gaussian Process (pGP). Both methods were applied to a hydraulic case corresponding to a spatially distributed open-channel steady flow along the Garonne River depending on the upstream discharge $Q$ and on the Strickler friction coefficient $K_{s_{3}}$, with the long-term objective to determine which surrogate strategy could be useful in the framework of ensemble-based data assimilation. It is important to show how surrogate models could be used to estimate some statistical quantities in a cost-effective way. This is useful for sensitivity analysis studies to evaluate the impact of physical parameters and external forcing on the river state. This is also useful for data assimilation to estimate correlation matrix and PDF related to the spatially distributed river state for the Ensemble Kalman Filter (EnKF) and the Particle Filter (PF), respectively.

We carried out a convergence study based on the following metrics: water level statistical moments, correlation matrix, PDF as well as Sobol’ indices representing the contribution of the upstream discharge and the Strickler friction coefficient on the water level variance. The accuracy of the PC and pGP surrogates were measured by their ability to retrieve the reference metrics obtained with a converged Monte Carlo (MC) random sample (including $\numprint{100,000}$ MASCARET simulations). An in-depth comparison of these metrics was done using the same computational budget: 121 MASCARET snapshots. The sensitivity to the number of snapshots was carried out to ensure that 121 simulations were enough to make this comparison valuable.

Results showed that both surrogate models can be used in place of the MASCARET forward model for uncertainty propagation without loss of accuracy. None of the two surrogate models clearly outperforms the other. In both cases, the PC and pGP surrogate models are able to correctly retrieve all physical information. The PC strategy seems to be more precise to compute the spatially distributed correlations as well as the Sobol’ indices, with the advantage that these indices do not need any further MASCARET evaluation as they are analytically computed from the PC expansion. Still, the multimodal water level PDF at Marmande (which is an important observation station along the Garonne River in operational context) was better captured by the pGP strategy that requires an additional sampling of the surrogate. Indeed, even increasing the number of snapshots to 256 and beyond was not enough to retrieve the second mode of the PDF using the PC surrogate model, while this was already captured with 121 snapshots by the pGP surrogate model. Still, it should be mentioned that the PC model better positions the first mode than the pGP model when using the MC approach as reference. Last but not least, the PC strategy requires some insight about the uncertain inputs of the system. We may not have access to this information in practice, leading potentially to a poor robustness of the PC surrogate. The quantity of interest may also feature non-linearity and exhibit extrema, which could be difficult to account for using quadrature points that could miss some physics. Alternative (for instance sparse) projection strategies could be investigated in the future to overcome these limitations.

Conclusions for the present test case highlight the validity of both quadrature-based PC and POD-based GP surrogate strategies for SWE in permanent flow for a small dimensional problem (the size of the uncertain space is $d=2$). The ranking between PC and pGP approaches will need to be further investigated when moving to open-channel unsteady flow modelling used for instance in the context of operational flood forecasting. The first challenge lies in the extension of the PC and pGP surrogates to larger uncertain dimension $d$, especially to address parameters that vary in space or over time or both such as the bathymetry spatial field and the time series of the upstream discharge. This may require the evaluation of more advanced strategies to reduce the size of the basis, the uncertain dimension (for instance through the Karhunen-Loève transformation) and the number of snapshots $N$ in the training set. For instance, the quadrature method used to build the PC surrogate is known to suffer from the curse-of-dimensionality. This will have to be revisited for a larger size $d$ of the uncertain space. The second challenge lies in the validity of the surrogate model over successive data assimilation cycles. It may be necessary to adjust the coefficients of the surrogates to track the changes in the river state behaviour over time. Those are crucial steps for the complementary use of surrogate models within data assimilation algorithms.