Stochastic Galerkin finite volume shallow flow model: well-balanced treatment over uncertain topography

TL;DR

This paper develops a stochastic Galerkin finite volume model for shallow water flows over uncertain topography, ensuring well-balanced solutions and significantly reducing computational cost compared to Monte Carlo methods.

Contribution

It introduces a low-order Wiener-Hermite Polynomial Chaos expansion with a finite volume approach that guarantees well-balancing and handles non-Gaussian uncertainties.

Findings

Model preserves lake-at-rest condition analytically and numerically.

Achieves comparable probability distributions to Monte Carlo with 100 times less computation.

Successfully models flows with discontinuous and non-Gaussian topography uncertainties.

Abstract

Stochastic Galerkin methods can quantify uncertainty at a fraction of the computational expense of conventional Monte Carlo techniques, but such methods have rarely been studied for modelling shallow water flows. Existing stochastic shallow flow models are not well-balanced and their assessment has been limited to stochastic flows with smooth probability distributions. This paper addresses these limitations by formulating a one-dimensional stochastic Galerkin shallow flow model using a low-order Wiener-Hermite Polynomial Chaos expansion with a finite volume Godunov-type approach, incorporating the surface gradient method to guarantee well-balancing. Preservation of a lake-at-rest over uncertain topography is verified analytically and numerically. The model is also assessed using flows with discontinuous and highly non-Gaussian probability distributions. Prescribing constant inflow over uncertain topography, the model converges on a steady-state flow that is subcritical or transcritical depending on the topography elevation. Using only four Wiener-Hermite basis functions, the model produces probability distributions comparable to those from a Monte Carlo reference simulation with 2000 iterations, while executing about 100 times faster. Accompanying model software and simulation data is openly available online.

Equations

$$\frac{\partial \mathbf{U}(x,t)}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{U}(x,t))}{\partial x}=\mathbf{S}(\mathbf{U}(x,t),z(x))$$

$${\mathbf{U}}_i^{(n+1)}={\mathbf{U}}_i^{(n)}-\Delta t\left(\frac{{\tilde{\mathbf{F}}}_{i+1/2}^{(n)}-{\tilde{\mathbf{F}}}_{i-1/2}^{(n)}}{\Delta x}-{\mathbf{S}}_i^{(n)}\right)$$

$${\mathbf{S}}_i^{(n)}={\left[0,-g\left(\frac{h_{i-1/2}^{+,\star }+h_{i+1/2}^{-,\star }}{2}\right)\left(\frac{z_{i+1/2}^{\star }-z_{i-1/2}^{\star }}{\Delta x}\right)\right]}^{⊺}$$

$$A(\theta )\approx \sum _{p=0}^P{A}_p{\Phi }_p(\xi (\theta ))$$

$${\Phi }_p(\xi )={\left(-1\right)}^p\exp \left(\frac{{\xi }^2}{2}\right)\frac{{\mathrm{d}}^p}{\mathrm{d}{\xi }^p}\exp \left(-\frac{{\xi }^2}{2}\right)$$

$$⟨{\Phi }_p{\Phi }_s⟩=⟨{\Phi }_p^2⟩{\delta }_{ps}$$

$$⟨\alpha (\xi )⟩={\int }_{-\infty }^{\infty }\alpha (\xi )W(\xi )\mathrm{d}\xi$$

$$W(\xi )=\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{{\xi }^2}{2}\right)$$

$$⟨{\Phi }_p{\Phi }_s⟩=⟨{\Phi }_s{\Phi }_p⟩$$

$$⟨A(\theta ){\Phi }_l⟩=A_l⟨{\Phi }_l^2⟩$$

$$⟨\left(A(\theta )+B(\theta )\right){\Phi }_l⟩=⟨A(\theta ){\Phi }_l⟩+⟨B(\theta ){\Phi }_l⟩$$

$${\mu }_m[A]={\int }_{-\infty }^{\infty }{\left(A-\beta \right)}^{m}W(\xi )\mathrm{d}\xi =⟨{\left(A-\beta \right)}^m⟩$$

$${\mu }_1[A]=A_0$$

$${\mu }_2[A]$$

$$f_A(a)=\sum _{j=1}^J{\left|\sum _{p=0}^P{A}_p\frac{\mathrm{d}{\Phi }_p}{\mathrm{d}\xi }({\hat{\xi }}_j)\right|}^{-1}W({\hat{\xi }}_j)$$

$$a-\sum _{p=0}^P{A}_p{\Phi }_p(\xi )=0$$

$${\mathbf{U}}_i^{(n)}(\theta )\approx \sum _{p=0}^P{\mathbf{U}}_{i,p}^{(n)}{\Phi }_p(\xi (\theta ))\text{,}{z}_i(\theta )\approx \sum _{p=0}^P{z}_{i,p}{\Phi }_p(\xi (\theta ))$$

$$\sum _{p=0}^P{z}_{i+1/2,p}^{\star }{\Phi }_p=\frac{1}{2}\left(\sum _{p=0}^P{z}_{i,p}{\Phi }_p+\sum _{p=0}^P{z}_{i+1,p}{\Phi }_p\right)$$

$$h_{i,p}^{(n)}={\eta }_{i,p}^{(n)}-z_{i,p}$$

$${\tilde{\mathbf{F}}}_{i+1/2}^{(n)}=\tilde{\mathbf{F}}\left(\sum _{p=0}^P{\mathbf{U}}_{i+1/2,p}^{-,\star }{\Phi }_p,\sum _{p=0}^P{\mathbf{U}}_{i+1/2,p}^{+,\star }{\Phi }_p\right)$$

$${\mathbf{S}}_i^{(n)}={\left[0,-\frac{g}{\Delta x}\left\{\sum _{p=0}^P\left(\frac{h_{i-1/2,p}^{+,\star }+h_{i+1/2,p}^{-,\star }}{2}\right){\Phi }_p\right\}\left\{\sum _{s=0}^P\left(z_{i+1/2,s}^{\star }-z_{i-1/2,s}^{\star }\right){\Phi }_s\right\}\right]}^{⊺}$$

$${\mathbf{U}}_{i,l}^{(n+1)}={\mathbf{U}}_{i,l}^{(n)}-\frac{\Delta t}{⟨{\Phi }_l^2⟩}\left(\frac{⟨{\tilde{\mathbf{F}}}_{i+1/2}^{(n)}{\Phi }_l⟩-⟨{\tilde{\mathbf{F}}}_{i-1/2}^{(n)}{\Phi }_l⟩}{\Delta x}-⟨{\mathbf{S}}_i^{(n)}{\Phi }_l⟩\right)$$

$$⟨{\tilde{\mathbf{F}}}_{i+1/2}^{(n)}{\Phi }_l⟩\approx \sum _{j=1}^{P+1}{w}_j\tilde{\mathbf{F}}\left(\sum _{p=0}^P{\mathbf{U}}_{i+1/2,p}^{-,\star }{\Phi }_p({\xi }_j),\sum _{p=0}^P{\mathbf{U}}_{i+1/2,p}^{+,\star }{\Phi }_p({\xi }_j)\right){\Phi }_l({\xi }_j)W({\xi }_j)$$

$$⟨{\mathbf{S}}_i^{(n)}{\Phi }_l⟩$$

$$\frac{⟨F_{i+1/2}^{(n)}{\Phi }_l⟩-⟨F_{i-1/2}^{(n)}{\Phi }_l⟩}{\Delta x}-⟨S_i^{(n)}{\Phi }_l⟩=0$$

$$F_{i+1/2}^{(n)}=\frac{g}{2}\left\{\sum _{p=0}^P{h}_{i+1/2,p}^{\star }{\Phi }_p\right\}\left\{\sum _{s=0}^P{h}_{i+1/2,s}^{\star }{\Phi }_s\right\}$$

$$⟨\frac{F_{i+1/2}^{(n)}-F_{i-1/2}^{(n)}}{\Delta x}{\Phi }_l⟩=⟨S_i^{(n)}{\Phi }_l⟩$$

$$⟨\frac{F_{i+1/2}^{(n)}-F_{i-1/2}^{(n)}}{\Delta x}{\Phi }_l⟩$$

$$h_{i-1/2,p}^{\star }{h}_{i+1/2,s}^{\star }-h_{i-1/2,p}^{\star }{h}_{i+1/2,s}^{\star })⟨{\Phi }_p{\Phi }_s{\Phi }_l⟩$$

$$⟨\frac{F_{i+1/2}^{(n)}-F_{i-1/2}^{(n)}}{\Delta x}{\Phi }_l⟩$$

$$-\sum _{p=0}^P\sum _{s=0}^P{h}_{i-1/2,p}^{\star }{h}_{i+1/2,s}^{\star }⟨{\Phi }_p{\Phi }_s{\Phi }_l⟩]$$

$$⟨\frac{F_{i+1/2}^{(n)}-F_{i-1/2}^{(n)}}{\Delta x}{\Phi }_l⟩$$

Full text

1 Introduction

Shallow water flows can be highly sensitive to uncertainties in input data: uncertainties in rating curve extrapolation, the choice of friction coefficient, digital elevation model errors and river survey measurement errors can substantially alter simulated distributions of flood extent, water depth and flow velocity (Bates et al., 2014; Kim and Sanders, 2016; Jung and Merwade, 2012; Montanari and Di Baldassarre, 2013). Conventional uncertainty quantification methods, including the Generalized Likelihood Uncertainty Estimation (GLUE) method (Beven and Binley, 1992), are based on Monte Carlo sampling with randomised inputs, making conventional methods computationally expensive and slow to converge on flow statistics. As a result, Monte Carlo simulations of flood events are severely constrained by available computing resources, and many sources of uncertainty must be neglected to make probabilistic simulations feasible (Neal et al., 2013). Stochastic Galerkin methods offer a computationally efficient alternative: repeated sampling is eliminated, and convergence is typically far more rapid than Monte Carlo methods (Xiu, 2009; Ge et al., 2008).

Polynomial Chaos methods work by introducing to the governing deterministic equations an additional, stochastic dimension, and approximating it with a so-called Polynomial Chaos expansion. Polynomial Chaos methods can be classified as non-intrusive or intrusive. Non-intrusive methods repeatedly sample a deterministic model with different input values, using the outputs to construct a stochastic solution. While non-intrusive polynomial chaos methods involve repeated sampling, they require far fewer samples than Monte Carlo methods (Ge et al., 2008; Ge and Cheung, 2011). Intrusive methods reformulate a deterministic model to produce a stochastic formulation that must be implemented as a new stochastic solver. The stochastic Galerkin method is an intrusive method which makes a Galerkin projection in stochastic space to produce a decoupled system of equations that are solved in a single model run to directly evolve the coefficients in the Polynomial Chaos expansion. Conventional stochastic Galerkin methods span the entire stochastic space using a single element with a basis chosen from the Askey scheme of orthogonal polynomials (Xiu and Karniadakis, 2002). When the stochastic flow is sufficiently smooth and the probability distributions are well-represented by the Polynomial Chaos basis, then Polynomial Chaos methods can achieve exponential convergence in the stochastic dimension (Xiu and Karniadakis, 2003). Sattar and El-Beltagy (2017) applied a conventional stochastic Galerkin method to the one-dimensional water hammer equations, and their model was able to converge on solutions with smooth, Gaussian probability distributions using a first-order Polynomial Chaos basis with just two basis functions.

Conventional stochastic Galerkin methods have two main shortcomings. First, strongly nonlinear flows produce steep gradients or discontinuities in stochastic space that are often poorly represented by global Polynomial Chaos bases using a single element (Pettersson et al., 2014). Second, even when exponential convergence is achieved early in the simulation, solution accuracy can deteriorate severely over long simulations (Gerritsma et al., 2010). These shortcomings have motivated a wealth of research into better alternatives, including time-dependent basis functions (Gerritsma et al., 2010), multiresolution wavelet decompositions (Le Maître et al., 2004), multi-element discretisations and adaptive meshing of stochastic space (Wan and Karniadakis, 2006; Tryoen et al., 2010; Pettersson et al., 2014; Li and Stinis, 2015).

While stochastic Galerkin methods have received a great deal of attention, few have applied the method to the shallow water equations. Ge et al. (2008) made the first step in this direction, formulating a one-dimensional stochastic Galerkin shallow water model to study solitary wave propagation over uncertain topography. Their stochastic Galerkin model was 50 times more efficient than Monte Carlo with similar accuracy, demonstrating that stochastic Galerkin shallow water models offer a viable alternative to Monte Carlo methods. However, their numerical tests were limited to flows with smooth, unimodal probability distributions that are more easily represented by stochastic Galerkin methods. Furthermore, their stochastic Galerkin formulation used a centred difference approximation of the bed slope source term that is not well-balanced.

A well-balanced stochastic Galerkin model was formulated by Jin et al. (2016) and applied to scalar nonlinear equations including the inviscid Burgers’ equation. The model inherited the well-balancing property from the underlying deterministic model that used the ‘interface method’ by Jin (2001). Jin et al. (2016) note that the intrusive stochastic Galerkin method has a particular advantage over non-intrusive stochastic collocation: a stochastic collocation method can only guarantee well-balancing at the collocation nodes where the deterministic well-balanced model is sampled. In contrast, the stochastic Galerkin model guarantees well-balancing across the entire stochastic space, irrespective of the order of Polynomial Chaos basis.

When applied to the shallow water equations, the well-balanced interface method is equivalent to the surface gradient method by Zhou et al. (2001). A stochastic reformulation of the surface gradient method is straightforward because it is a linear method involving only linear combinations of the discrete flow variables and topography. Popular alternatives to the surface gradient method include well-balanced approaches that are also depth-positivity preserving (Audusse et al., 2004; Liang and Marche, 2009). Such approaches use nonlinear $\max(\cdot,\cdot)$ operators which greatly complicate their stochastic Galerkin reformulation for realistic applications.

This paper presents a well-balanced stochastic Galerkin model of the one-dimensional shallow water equations for simulating probabilistic flows that account for measurement errors in flood plain topography or open channel bathymetry. This one-dimensional stochastic model represents the first step towards a two-dimensional probabilistic flood model, and it is also relevant to one-dimensional open channel flows. Being based on a Godunov-type shock-capturing method, the model is capable of simulating subcritical, supercritical and transcritical flows. The model uses a stochastic reformulation of the surface gradient method (Zhou et al., 2001), and the well-balancedness of this reformulation is studied theoretically and verified numerically for an idealised, motionless lake-at-rest over an uncertain and irregular bed. Next, a steady-state test is proposed that provokes a strongly nonlinear flow response over an uncertain bed, resulting in probability distributions that are highly non-Gaussian and discontinuous. Probability distributions produced by the stochastic Galerkin model are validated against a Monte Carlo reference simulation, and the relative computational performance is discussed. Finally, the stochastic Galerkin model is verified for a more realistic steady-state flow over a highly irregular bed with uncertainties that characterise measurement errors in bed elevation. Numerical simulation data (Shaw and Kesserwani, 2019a) and a Python 3 implementation of the stochastic Galerkin shallow flow model (Shaw and Kesserwani, 2019b) are available to download from Zenodo. Instructions for running the models and interpreting the data are provided in the appendix.

2 Deterministic and stochastic shallow flow models

In this section, a certain deterministic numerical solver of the one dimensional (1D) shallow water equations is outlined in the framework of a finite volume Godunov-type method. The selected deterministic solver relies on the surface gradient method (Zhou et al., 2001) to both ensure a well-balanced topography integration and extendibility of the well-balanced property into the stochastic Galerkin case. Accordingly, a stochastic Galerkin reformulation is devised that is theoretically well-balanced with uncertain topography under a lake-at-rest hypothesis.

The mathematical model of the shallow water equations represent mass and momentum conservation principals, and is used in the following conservative form when solving it within a finite volume Godunov-type framework (Toro and Garcia-Navarro, 2007):

[TABLE]

where $\mathbf{U}=\left[h,q\right]^{\intercal}$ is the flow vector including water depth $h$ ($\mathrm{L}$) and unit-width discharge $q=hv$ ($\mathrm{L}^{2}/\mathrm{T}$) in which $v$ represents the depth-averaged velocity ($\mathrm{L}/\mathrm{T}$), $\mathbf{F}=\left[q,q^{2}/h+gh^{2}/2\right]^{\intercal}$ is the flux vector in which $g$ represents the gravitational constant and $\mathbf{S}=\left[0,-gh\>\mathrm{d}z/\mathrm{d}x\right]^{\intercal}$ is the source term vector in which the gradient of the topography $z(x)$ is involved. Equation (1) represents hydraulic flow in an idealised, frictionless channel with a rectangular cross-section of unit width.

2.1 Deterministic model

On a uniform 1D mesh with $M$ elements each of size $\Delta x$, the first-order finite volume method leads to the following discrete element-wise formulation of the shallow water equations given by equation (1):

[TABLE]

in which $\mathbf{U}_{i}^{(n)}=\left[h_{i}^{(n)},q_{i}^{(n)}\right]^{\intercal}$ is a piecewise-constant discretisation of the flow vector at element $i$ and time level $(n)$, and $\mathbf{\tilde{F}}_{i+1/2}^{(n)}$ is a numerical flux function for linking nonlinear discontinuities associated with the flow vector data at interface $i+1/2$ located between element $i$ and $i+1$. Namely, $\mathbf{\tilde{F}}_{i+1/2}^{(n)}=\mathbf{\tilde{F}}(\mathbf{U}_{i+1/2}^{-},\mathbf{U}_{i+1/2}^{+})$ where $\mathbf{U}_{i+1/2}^{-}$ is the limit of the solution from the side of element $i$, and $\mathbf{U}_{i+1/2}^{+}$ is the limit of the solution from the side of element $i+1$. Within the scope of this work involving a first-order accurate solver, these limits become $\mathbf{U}_{i+1/2}^{-}=\mathbf{U}_{i}$ and $\mathbf{U}_{i+1/2}^{+}=\mathbf{U}_{i+1}$ that are used in a numerical flux function based on the Roe approximate Riemann solver (Roe and Pike, 1984). Consequently, the deterministic model is able to capture the occurrence of shocks, and simulate subcritical, supercritical and transcritical flows.

The surface gradient method essentially reconstructs an averaged topography at interface $i+1/2$ that is shared by both elements $i$ and $i+1$, as $z^{\star}_{i+1/2}=(z_{i}+z_{i+1})/2$. From the reconstructed topography $z^{\star}_{i+1/2}$, consistent flow variable limits are accordingly reconstructed based on the actual free-surface elevation data, i.e. $\eta_{i+1/2}^{-}=h_{i}^{(n)}+z_{i}$ and $\eta_{i+1/2}^{+}=h_{i+1}^{(n)}+z_{i+1}$, and velocity data, i.e. $v_{i+1/2}^{-}=q_{i}^{(n)}/h_{i}^{(n)}$ and $v_{i+1/2}^{+}=q_{i+1}^{(n)}/h_{i+1}^{(n)}$, as: $h^{K,\star}_{i+1/2}=\eta_{i+1/2}^{K}-z^{\star}_{i+1/2}$ and $q^{K,\star}_{i+1/2}=h^{K,\star}_{i+1/2}v_{i+1/2}^{K}$ (where $K=+\text{ or }-$). For clarity of presentation, the time level denoted by superscript $(n)$ is omitted from all reconstructed variables. These reconstructions form new Riemann states $\mathbf{U}^{K,\star}_{i+1/2}=\left[h^{K,\star}_{i+1/2},q^{K,\star}_{i+1/2}\right]^{\intercal}$ for use to evaluate $\mathbf{\tilde{F}}_{i+1/2}^{(n)}$. By analogy, new Riemann limits $\mathbf{U}^{K,\star}_{i-1/2}=\left[h^{K,\star}_{i-1/2},q^{K,\star}_{i-1/2}\right]^{\intercal}$ at $i-1/2$ are produced for use to evaluate $\mathbf{\tilde{F}}_{i-1/2}^{(n)}$. From the reconstructed limits, a well-balanced discretisation of the source term vector can be produced:

[TABLE]

The well-balanced deterministic model presented in equations (2) and (3) is used for Monte Carlo simulations, and the deterministic model is also the starting point for a stochastic Galerkin reformulation.

2.2 Fundamental properties of the Polynomial Chaos basis

Before presenting the stochastic Galerkin reformulation, it is necessary to consider a single random variable $A(\theta)$ that maps from the random event $\theta$ to an arbitrary probability distribution with finite variance. This random variable can be approximated by a Wiener-Hermite Polynomial Chaos expansion (Xiu and Karniadakis, 2002). The expansion is based on a standard Gaussian random variable $\xi(\theta)\in[-\infty,+\infty]$ having zero mean and unit variance. The random variable of interest, $A(\theta)$, is then approximated as

[TABLE]

where $\mathbf{A}=\left[A_{0},\ldots,A_{P}\right]^{\intercal}$ are the expansion coefficients and $\bm{\Phi}=\left[\Phi_{0},\ldots,\Phi_{P}\right]^{\intercal}$ is the probabilists’ Hermite polynomial basis having basis function $\Phi_{p}$ of degree $p$,

[TABLE]

where $\Phi_{0}=1,\Phi_{1}=\xi,\Phi_{2}=\xi^{2}-1,\Phi_{3}=\xi^{3}-3\xi$ and so on. As the basis order $P$ is increased, the Wiener-Hermite Polynomial Chaos approximation converges on the true random variable $A(\theta)$ (Xiu and Karniadakis, 2002).

Basis orthogonality and commutativity

The Wiener-Hermite basis $\bm{\Phi}$ is orthogonal such that

[TABLE]

where $\left\langle\cdot\right\rangle$ is the ensemble average operator and $\delta_{ps}$ is the Kronecker delta that is equal to one when $p=s$ and zero otherwise. The ensemble average operator is defined as the weighted integral over the standard Gaussian random variable $\xi$:

[TABLE]

where $\alpha(\xi)$ is an expression involving any combination of random variables or basis functions, and the weighting function $W(\xi)$ is the standard Gaussian probability density function

[TABLE]

This weighting function ensures that, when $\alpha$ is independent of $\xi$, the ensemble average $\left\langle\alpha\right\rangle=\alpha$. Finally, the ensemble average of a product of basis functions is commutative such that

[TABLE]

The commutative property is needed later when verifying the well-balanced property with uncertain topography.

Stochastic Galerkin projection of a random variable

Given the random variable $A(\theta)=\sum_{p=0}^{P}A_{p}\Phi_{p}(\xi(\theta))$, its Galerkin projection onto a basis function $\Phi_{l}$ with $l=0,\ldots,P$ is achieved using the ensemble average operator such that, due to orthogonality,

[TABLE]

where $A_{l}$ is the $l$th order expansion coefficient. Also note that the Galerkin projection of a basis function $\Phi_{l}$ and two random variables, $A(\theta)$ and $B(\theta)$, is distributive:

[TABLE]

Mean, variance and high-order moments

The mean, variance and high-order moments can be calculated for $A(\theta)$. The $m$th moment $\mu_{m}[A]$ is defined as

[TABLE]

where $\beta=0$ when $m=1$ and $\beta=\mu_{1}[A]$ for higher-order moments. Therefore, the mean $\mu_{1}[A]=\left\langle A\right\rangle=\sum_{p=0}^{P}A_{p}\left\langle\Phi_{p}\right\rangle$. Since $\left\langle\Phi_{0}\right\rangle=1$ and $\left\langle\Phi_{p}\right\rangle=0$ for $p>0$ then

[TABLE]

The shorthand notation for the mean of $A$ is $\overline{A}$, also known as the expected value, $\operatorname{\mathbb{E}}\left[A\right]$.

The variance $\mu_{2}[X]$ can be derived using the fact that $\operatorname{\mathbb{E}}\left[\left(A-\operatorname{\mathbb{E}}[A]\right)^{2}\right]=\operatorname{\mathbb{E}}[A^{2}]-\operatorname{\mathbb{E}}^{2}[A]$, hence $\mu_{2}[A]=\left(\sum_{p=0}^{P}A_{p}^{2}\left\langle\Phi_{p}^{2}\right\rangle\right)-A_{0}^{2}\left\langle\Phi_{0}\right\rangle^{2}$. Since $\left\langle\Phi_{0}\right\rangle^{2}=\left\langle\Phi_{0}^{2}\right\rangle$ then

[TABLE]

The shorthand notation for the variance of $A$ is $\sigma^{2}_{A}$ and the standard deviation of $A$ is $\sigma_{A}$.

Reconstructing the probability density function

The probability density function $f_{A}(a)$ of a random variable $A$ is,

[TABLE]

which can be calculated numerically for a specific realisation $a$. Hence, the probability density function is computed by evaluating equation (15) for a range of outcomes.

2.3 Stochastic Galerkin reformulation of the deterministic model

The solution of the stochastic 1D shallow water equations is now random because it depends on uncertain initial conditions, uncertain boundary conditions and uncertain topography. Hence, the stochastic 1D shallow water equations depend not only upon space $x$ and time $t$, but additionally upon the random event $\theta$. The stochastic flow vector $\mathbf{U}(x,t,\theta)$ becomes a general stochastic process having arbitrary probability distributions that vary in space and time. Similarly, the stochastic topography $z(x,\theta)$ has arbitrary probability distributions that vary in space.

The stochastic Galerkin reformulation of the deterministic model involves three steps to (i) replace the deterministic variables, $\mathbf{U}_{i}^{(n)}$ and $z_{i}$, with random variables $\mathbf{U}_{i}^{(n)}(\theta)$ and $z_{i}(\theta)$, (ii) rewrite the deterministic formulation using these random variables, and (iii) make a stochastic Galerkin projection onto the Wiener-Hermite basis.

Replacing deterministic variables with random variables

For all elements $i=1,\ldots,M$ across all time levels, every deterministic flow variable $\mathbf{U}_{i}^{(n)}=\left[h_{i}^{(n)},q_{i}^{(n)}\right]^{\intercal}$ and deterministic topography variable $z_{i}$ becomes a random variable approximated by a Wiener-Hermite Polynomial Chaos expansion:

[TABLE]

where $\mathbf{U}_{i,p}^{(n)}=\left[h_{i,p}^{(n)},q_{i,p}^{(n)}\right]^{\intercal}$ and $z_{i,p}$ are the $p$th order expansion coefficients over element $i$ at time level $n$.

The reconstructed topography and reconstructed limits become functions of random variables. The reconstructed topography at interface $i+1/2$ becomes

[TABLE]

and so $z^{\star}_{i+1/2,p}=(z_{i,p}+z_{i+1,p})/2$ due to basis orthogonality. The reconstructed limits $\mathbf{U}^{K,\star}_{i+1/2,p}=\left[h^{K,\star}_{i+1/2,p},q^{K,\star}_{i+1/2,p}\right]^{\intercal}$ (where $K=+\text{ or }-$) are calculated in a similar fashion.

Random variables that are functions of other random variables can be calculated in the same way. In particular, water depth can be expressed as a function of free-surface elevation and topography such that, due to basis orthogonality,

[TABLE]

which is used later for specifying initial conditions.

Rewriting the deterministic formulation using random variables

The deterministic finite volume formulation given by equation (2) is rewritten in terms of the random variables in equation (16). As a result, the numerical fluxes $\mathbf{\tilde{F}}_{i+1/2}^{(n)}$ and $\mathbf{\tilde{F}}_{i-1/2}^{(n)}$ and source term vector $\mathbf{S}_{i}^{(n)}$ become functions of random variables. The numerical flux $\mathbf{\tilde{F}}_{i+1/2}^{(n)}$ becomes

[TABLE]

and similarly for $\mathbf{\tilde{F}}_{i-1/2}^{(n)}$. The source term vector $\mathbf{S}_{i}^{(n)}$ in equation (3) becomes

[TABLE]

Equation (20) involves the product of two expressions, each delimited by braces. Since both expressions include Wiener-Hermite expansions then different indices, $p$ and $s$, are needed for the expansion coefficients in each expression.

Stochastic Galerkin projection

Due to the orthogonal property (equation 10) and distributive property (equation 11) of the Wiener-Hermite basis, a Galerkin projection of equation (2) onto the basis functions $\Phi_{l},l=0,\ldots,P$ produces $P+1$ decoupled equations:

[TABLE]

Equation (21) involves ensemble averages of numerical fluxes, $\left\langle\mathbf{\tilde{F}}_{i+1/2}^{(n)}\Phi_{l}\right\rangle$ and $\left\langle\mathbf{\tilde{F}}_{i-1/2}^{(n)}\Phi_{l}\right\rangle$, and an ensemble average of the source term vector, $\left\langle\mathbf{S}_{i}^{(n)}\Phi_{l}\right\rangle$. There is no straightforward method for calculating an ensemble average of the numerical flux because it is nonlinear. Instead, the integral in equation (7) is approximated by Gauss-Hermite quadrature,

[TABLE]

where $w_{j}$ are the quadrature weights and $\xi_{j}$ are the quadrature points. The ensemble average $\left\langle\mathbf{\tilde{F}}_{i-1/2}^{(n)}\Phi_{l}\right\rangle$ is calculated in the same way.

Unlike the nonlinear numerical flux, the ensemble average of the source term vector $\left\langle\mathbf{S}_{i}^{(n)}\Phi_{l}\right\rangle$ is linear and can be derived directly from equation (20):

[TABLE]

The ensemble averages $\left\langle\Phi_{p}\Phi_{s}\Phi_{l}\right\rangle$ in equation (23) and $\left\langle\Phi_{l}^{2}\right\rangle$ in equation (21) can be calculated analytically or exactly by Gauss-Hermite quadrature. Since these calculations do not depend on the solution then they can be precomputed once and stored.

2.4 Well-balanced Property with Uncertain Topography

In this section it is shown that the stochastic Galerkin model satisfies the well-balanced property for a lake-at-rest with uncertain topography. A lake at rest is an idealised scenario where the free-surface is entirely flat and motionless, much like still water in a lake. Assuming a lake-at-rest for the stochastic model, the mean free-surface elevation is constant and the mean discharge is zero, and there is no uncertainty in the free-surface elevation or discharge. The bed elevation can have any spatial profile and stochastic profile.

Since the free-surface elevation is constant and the averaged bed elevation is continuous at interfaces then $\eta_{i+1/2,p}^{-}-z^{\star}_{i+1/2,p}=\eta_{i+1/2,p}^{+}-z^{\star}_{i+1/2,p}$, hence $h^{-,\star}_{i+1/2,p}=h^{+,\star}_{i+1/2,p}$ for all $i=0,\ldots,M$ and $p=0,\ldots,P$. Shorthand notation is introduced such that $\eta_{p}$ is the $p$th expansion coefficient of the spatially uniform free-surface elevation, and $h^{\star}_{i+1/2,p}=h^{-,\star}_{i+1/2,p}=h^{+,\star}_{i+1/2,p}$. Since $h^{\star}_{i+1/2,p}$ is continuous at interfaces and the discharge is zero then the numerical flux is equal to the physical flux.

In order to preserve a lake-at-rest solution in element $i$ then the ensemble average of the flux gradient must balance the ensemble average of the source term vector in equation (21). The balance of mass and momentum components can be considered separately. The mass continuity equation balances because the discharge is zero. For the momentum equation, the following equality must hold for all $l=0,\ldots,P$:

[TABLE]

where the ensemble average of the bed slope source term $\left\langle S_{i}^{(n)}\Phi_{l}\right\rangle$ is given by equation (23) and the momentum flux $F_{i+1/2}^{(n)}$ is

[TABLE]

and similarly for $F_{i-1/2}^{(n)}$. Due to the distributive property of the basis (equation 11), equation (24) can be rewritten as

[TABLE]

The ensemble average of the momentum flux gradient is

[TABLE]

Equation (32) is equal to the ensemble average of the bed slope source term $\left\langle S_{i}^{(n)}\Phi_{l}\right\rangle$ given in equation (23), hence discrete balance is preserved.

3 Idealised numerical tests

Three numerical tests are performed to validate the well-balanced stochastic Galerkin model. The first test simulates a lake-at-rest over idealised, uncertain topography to verify that the stochastic Galerkin model preserves well-balancing numerically. The second test simulates flow over a hump with an uncertain elevation. This test is designed to challenge the stochastic Galerkin model at representing discontinuous, non-Gaussian probability distributions by generating a steady-state solution that may be subcritical or transcritical depending on the hump elevation. The third test verifies the robustness of the stochastic Galerkin model for a steady flow over a highly irregular and uncertain bed that is more representative of real-world river hydraulics.

3.1 Specification of Tests 1 and 2

For the first two tests, the 1D domain is [$-50\text{\,}\mathrm{m}$, $50\text{\,}\mathrm{m}$], tessellated by $M=100$ elements with no overlaps or gaps such that the mesh spacing is $\Delta x=$1\text{,}\mathrm{m}. The timestep is $\Delta t=$0.15\text{\,}\mathrm{s} resulting in a maximum Courant number of about $0.8$. By choosing a fixed timestep, simulations of a given test complete in the same number of timesteps irrespective of the model configuration, and error accumulation due to timestepping errors will be the same across all models.

Both tests include a topographic hump centred at $x=$0\text{,}\mathrm{m}$$ with a region of Gaussian uncertainty. Following a similar approach to Ge et al. (2008), there are two representations for the same uncertain topography. The first representation enables smooth topography profiles to be randomly generated in Monte Carlo iterations. The topography $z$ is defined as

[TABLE]

where the hump amplitude $r$ is a random variable with mean $\overline{r}=$0.6\text{,}\mathrm{m} and standard deviation $\sigma_{r}=$0.3\text{\,}\mathrm{m}, and the half-width is $\lambda=$10\text{,}\mathrm{m}$$. This topography profile is seen in figure 3. The second representation is used by the stochastic Galerkin model, with topography represented by expansion coefficients $z_{i,0},\ldots,z_{i,P}$. To be able to calculate the topography expansion coefficients, equation (33) must be expressed in terms of mean topography $\overline{z}(x)$ and topographic variance $\sigma_{z}^{2}(x)$ without involving the random variable $r$. The mean topography $\overline{z}(x)$ is simply

[TABLE]

The topographic variance is

[TABLE]

Equation (35) can be rewritten using Taylor series expansions of the two terms $\operatorname{\mathbb{E}}\left[z^{2}(x,r)\right]$ and $\operatorname{\mathbb{E}}^{2}\left[z(x,r)\right]$. To illustrate the approach, a Taylor series expansion of $\operatorname{\mathbb{E}}\left[z(x,r)\right]$ about $\overline{r}$ is

[TABLE]

Applying this approach to equation (35) gives:

[TABLE]

Equipped with analytic expressions for the mean topography $\overline{z}(x)$ (equation 34) and topographic variance $\sigma_{z}^{2}(x)$ (equation 39), now the topography expansion coefficients $z_{i,0},\ldots,z_{i,P}$ can be calculated. Since the topographic bump has a Gaussian probability distribution with $\mu_{1}[z(x)]=\overline{z}(x)$, $\mu_{2}[z(x)]=\sigma_{z}^{2}(x)$ and high-order moments $\mu_{m}[z(x)]=0$ for $m\geq 3$ then, using equations (12), (13) and (14), the topography expansion coefficients are

[TABLE]

where values are calculated at the centre point $x_{i}$ for all elements $i=1,\ldots,M$.

The initial water depth expansion coefficients $h_{i,0}^{(0)},\ldots,h_{i,P}^{(0)}$ can be calculated in terms of free-surface elevation and topography using equation (18). For both tests, the initial, spatially-uniform mean free-surface elevation is $1.5\text{\,}\mathrm{m}$ with no initial uncertainty such that $\eta_{i,0}^{(0)}=$1.5\text{,}\mathrm{m}$$ and $\eta_{i,p}^{(0)}=0$ with $p=1,\ldots,P$ and $i=1,\ldots,M$.

3.2 Test 1: Lake-at-Rest Over an Uncertain Bed

Numerical methods that are not well-balanced produce spurious waves in the vicinity of sloping topography, and these spurious waves are particularly evident for slow-moving flows with weak momentum fluxes. In the limit, the momentum flux is zero and the water is motionless, like the water surface on a lake at rest. Hence, the lake-at-rest test is ideally suited to verify the well-balanced property, since the analytic solution preserves the resting state forever, and any waves generated by a numerical model are entirely spurious (Bermúdez and Vázquez, 1994).

To present a challenging test, a rectangular obstacle is introduced to the right of the uncertain hump given by equation (33), so the bed elevation $z$ becomes

[TABLE]

Results of the well-balanced stochastic model are compared with those of a stochastic model having a centred difference approximation of the source term vector that does not exactly balance the numerical flux gradient. The centred difference model is the same as the well-balanced stochastic model except for two changes. First, numerical fluxes are calculated using the original, unmodified flow variables:

[TABLE]

Second, the ensemble average of the source term vector uses a centred difference approximation:

[TABLE]

The centred difference model and the well-balanced model are both configured with a Wiener-Hermite basis order $P=3$.

The simulated time is $100\text{\,}\mathrm{s}$ corresponding to about 670 timesteps, and the solutions for the centred difference and well-balanced stochastic Galerkin models are shown in figure 1. The lack of well-balancing is apparent using the centred difference model: grid-scale standing waves develop at the discontinuities either side of the rectangular hump (figure 1, 1), and a smooth standing wave also develops over the uncertain hump. These errors persist throughout the simulation. In contrast, the well-balanced stochastic Galerkin model preserves the initial resting state with discharges accurate to machine precision (figure 1). This numerical result confirms that the stochastic Galerkin model is well-balanced in theory and in practice.

The choice of the Wiener-Hermite basis introduces a particular limitation that imposes an upper bound on the basis order $P$, and constrains the minimum water depth that the stochastic Galerkin model can represent. This limitation arises because the hump amplitude has a Gaussian probability distribution so the tails of the distribution extend to $\pm\infty$, meaning that there is a non-zero probability that the water depth is negative. The stochastic Galerkin formulation presented here does not accommodate wetting-and-drying processes, and any negative water depth will crash the model. If the basis order $P$ is increased then the Gauss-Hermite quadrature points in equation (22) extend further into the tails of the probability distributions, leading to negative water depths being provided as input to the Riemann solver. Similarly, raising the topography, decreasing the initial water depth, or increasing the topographic uncertainty can all produce negative water depths in the stochastic Galerkin model. This behaviour has been verified experimentally by varying the model basis order, initial conditions and topography profile.

3.3 Test 2: Steady-State Critical Flow Over an Uncertain Bed

This test is designed to challenge the stochastic Galerkin model at representing highly non-Gaussian and discontinuous distributions of stochastic flow. The uncertain topography profile and inflow boundary condition are chosen in order to produce a nonlinear response such that the steady-state solution may be subcritical or transcritical depending on the bed elevation. Results of the stochastic Galerkin model are validated against a Monte Carlo simulation that serves as a reference solution.

3.3.1 Setup to Produce a Nonlinear Flow Response

The uncertain topography is the smooth hump given by equation (40). Subcritical boundary conditions are imposed such that the mean upstream discharge per unit-width is $1.65\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ and the mean downstream water depth is $1.5\text{\,}\mathrm{m}$, with no uncertainty for the upstream discharge and downstream water depth. Transmissive boundary conditions are used for the upstream water depth and downstream discharge. These boundary conditions are chosen so that the flow is exactly critical for the mean hump amplitude $\overline{r}=$0.6\text{,}\mathrm{m} at $x=$0\text{\,}\mathrm{m}. Since the hump amplitude is uncertain then the flow regime is also uncertain: if the hump amplitude is less than $\overline{r}$ then the flow remains subcritical; if the hump amplitude is greater than $\overline{r}$ then the flow regime becomes transcritical. In the transcritical regime, the flow upstream of the hump is subcritical, transitioning to supercritical flow over the hump. A hydraulic jump occurs on the downstream side of the hump, where the flow becomes subcritical once more.

To illustrate this change in flow regime, figure 2 shows four deterministic solutions using four different hump amplitudes. Solutions from the well-balanced deterministic model are obtained at $t=$500\text{,}\mathrm{s}$$ when the water has converged on a steady state. Convergence is measured by calculating the $L^{2}$ difference in mean water depth between the current and previous timesteps:

[TABLE]

By $t=$500\text{,}\mathrm{s} all four deterministic solutions have converged down to a convergence error of $10^{-4}$ $\mathrm{m}$. For a small hump with amplitude $\overline{r}-\sigma_{r}=$0.3\text{\,}\mathrm{m}, the flow remains subcritical. A linear increase in hump amplitude produces a strongly nonlinear response in the steady-state water profile, as seen in figure 2. Two nonlinear responses are evident in particular: first, the upstream boundary condition allows the upstream water depth to increase nonlinearly; second, a transcritical shock develops over the hump that increases in amplitude and moves further downstream with larger hump amplitudes. Downstream of the hump, the water depth is $1.5\text{\,}\mathrm{m}$ irrespective of the hump amplitude, with this profile having propagated upstream from the imposed downstream boundary.

3.3.2 Configuration of the Monte Carlo Reference Simulation

The Monte Carlo reference simulation is performed using iterations of the well-balanced deterministic model. It is necessary to perform a sufficient number of Monte Carlo iterations to ensure that the flow statistics are accurate, and more iterations are needed for more complex probability distributions. After each Monte Carlo iteration, the mean and standard deviation of water depth at $x=$1.5\text{,}\mathrm{m}$$ are measured, where the probability distribution is most complex. The statistics are compared with those of the previous iteration, and the Monte Carlo simulation terminates once the statistics change negligibly between iterations, when statistical convergence is achieved. Guided by these measurements, two thousand Monte Carlo iterations are needed to achieve statistical convergence for this test.

For each Monte Carlo iteration, the topography is randomly generated using a hump amplitude drawn from the Gaussian distribution given by $(\overline{r},\sigma_{r})$ and so the topography will always be smooth. If instead the topography was randomly generated using $(\overline{z}(x),\sigma_{z}(x))$ then randomisation would be different in every element and the topography would not be smooth, so many more iterations would be needed to sample the stochastic solution space. For the Monte Carlo iterations, the hump amplitude $r$ is constrained such that 0\text{\,}\mathrm{m}$\leq r\leq$1.4\text{\,}\mathrm{m} to avoid negative water depths.

3.3.3 Spatial Profiles of the Uncertain Free-Surface Elevation

In figure 3, spatial profiles of the uncertain free-surface elevation are obtained at $t=$500\text{,}\mathrm{s}$$ when the water depth profiles from the Monte Carlo and stochastic Galerkin simulations have converged down to $10^{-4}$ $\mathrm{m}$ as defined by equation (44). Using a Wiener-Hermite basis of order $P=3$, the stochastic Galerkin model accurately represents the Monte Carlo reference profiles of the mean and standard deviation of free-surface elevation. Upstream of the hump, the stochastic Galerkin model predicts a standard deviation that is slightly too small compared to the Monte Carlo reference solution. Small errors in the stochastic Galerkin free-surface elevation are also visible above the hump where the flow is most complex. Similar results are obtained using basis order $P=1$ or $P=2$, with stochastic Galerkin errors increasing slightly as the basis order is decreased (not shown).

3.3.4 Monte Carlo Histograms of Uncertain Free-Surface Elevation

The mean and standard deviation statistics are useful for summarising the spatial profile of uncertainty, but they are less meaningful for non-Gaussian probability distributions. In such cases, it is more meaningful to study the complete probability distributions, which is also particularly important for flood risk assessments that are concerned with extreme events that occur in the tails of the distributions (Ge and Cheung, 2011).

Probability distributions of the free-surface elevation are sampled at two points at $t=$500\text{,}\mathrm{s}: the first at $x=$-37.5\text{\,}\mathrm{m} and the second at $x=$1.5\text{,}\mathrm{m}$$, with these positions marked by dotted lines in figure 2 and figure 3. The first point is far upstream of the hump where the free-surface elevation is uncertain and spatially uniform. The second point is immediately downstream of the hump in the region where transcritical shocks can develop.

Figure 4 shows histograms from the Monte Carlo reference simulation that estimate the true probability densities at the two points. Stochastic Galerkin results which appear on the same figure are discussed later. For each of the two points, water depths from the 2000 Monte Carlo iterations are binned into intervals of $0.05\text{\,}\mathrm{m}$, and the magnitude of each bin represents the probability that the water depth is within the given interval. Since the histogram estimates the probability density then the total shaded area over all bins is equal to one.

The Monte Carlo histogram at $x=$-37.5\text{,}\mathrm{m} is shown in figure [4](#S3.F4) and is discussed first. For subcritical flows over small humps, the upstream water level remains at its initial height of $1.5\text{\,}\mathrm{m}$. For transcritical flows over larger humps, the upstream water level increases nonlinearly. Since the initial conditions and boundary conditions are chosen so that the mean flow is critical, then about 50% of the flows are subcritical, resulting in a large peak in the probability distribution at $\eta=$1.5\text{\,}\mathrm{m}. The other 50% of the flows are transcritical with elevated upstream water levels, resulting in a long tail in the distribution.

The stochastic flow response immediately downstream of the hump at $x=$1.5\text{,}\mathrm{m} is more complex. For subcritical flows, the steady-state free-surface elevation at this point will be only slightly lower than the initial free-surface elevation of $1.5\text{\,}\mathrm{m}$. For transcritical flows, the steady-state free-surface elevation may be much lower since, at $x=$1.5\text{\,}\mathrm{m}, the transcritical shock is close to its minimum depth. The subcritical and transcritical regimes appear as a bimodal distribution in the histogram in figure 4 with one peak around $\eta=$1.0\text{,}\mathrm{m} associated with transcritical flow and a second peak around $\eta=$1.5\text{\,}\mathrm{m} associated with subcritical flow.

3.3.5 Stochastic Galerkin Free-Surface Elevation Probability Densities

Overlayed on the Monte Carlo histograms, figure 4 also shows probability density functions obtained from stochastic Galerkin simulations. Three stochastic Galerkin simulations are performed using Wiener-Hermite bases of order $P=1$, $P=2$ and $P=3$. Free-surface elevation expansion coefficients are calculated by rearranging equation (18), from which probability density functions are reconstructed using equation (15). If a probability density function was in exact agreement with a Monte Carlo histogram then the line would pass through the top of every histogram bin, and so any deviation represents a numerical error associated with the stochastic Galerkin model.

Using basis order $P=1$, the stochastic Galerkin model can only represent Gaussian distributions, with its two expansion coefficients representing the mean and standard deviation. While figure 3 confirms that the mean and standard deviation statistics from the stochastic Galerkin model are accurate, the Monte Carlo histograms cannot be well-represented by Gaussian distributions.

Increasing the basis order to $P=2$, the probability distribution of upstream water levels is in good agreement with the Monte Carlo histogram (figure 4), though an error is noticeable at the discontinuity: the Monte Carlo histogram has a discontinuity at $\eta=$1.5\text{,}\mathrm{m} because the upstream water level can only rise above its initial height, so the probability that $\eta<$1.5\text{\,}\mathrm{m} is zero. The stochastic Galerkin model with basis order $P=2$ slightly underestimates this discontinuity, placing it at about $\eta=$1.45\text{,}\mathrm{m}. The distribution in free-surface elevation at $x=$1.5\text{\,}\mathrm{m} (figure 4) is not improved using basis order $P=2$, with the distribution remaining close to Gaussian.

Finally, the basis order is increased to $P=3$. The distribution of upstream water levels shifts slightly to the right in figure 4, corresponding to slightly higher water levels. As a result, the discontinuity is closer to the true value of $\eta=$1.5\text{,}\mathrm{m}, but this shift also produces slightly larger errors in the tail of the distribution. The most notable improvement is seen in the distribution at $x=$1.5\text{\,}\mathrm{m} (figure 4): the distribution from the stochastic Galerkin model now has a similar shape to the Monte Carlo histogram. The lower bound at $\eta=$1.0\text{,}\mathrm{m} is accurately represented, but the true upper bound around $\eta=$1.5\text{\,}\mathrm{m} is overestimated by about $0.1\text{\,}\mathrm{m}$. The probability density function has singularities at the lower and upper bounds that appear as spikes in the plot. These two singularities occur because the probability density function is the derivative of the cumulative density function, which is discontinuous and non-differentiable at these points.

While the stochastic Galerkin model with basis order $P=3$ adequately represented the true distribution bounds and distribution shapes, the probability density functions were not entirely accurate. Such inaccuracies are to be expected because low-order Wiener-Hermite bases cannot represent complex distributions. The basis order cannot be increased beyond $P=3$ because water depths are small near transcritical shocks, and the stochastic Galerkin model crashes for the same reason as discussed in the lake-at-rest test. Instead, results might be improved by choosing a more sophisticated method to discretise stochastic space. One candidate is the stochastic Galerkin multiwavelet approach by Pettersson et al. (2014), which is able to simulate stochastic gas dynamics with densities close to zero, analogous to very small water depths in shallow water flows.

3.3.6 Storage Requirements and Computation Time

Since the stochastic Galerkin model is necessarily more complex than the deterministic model, associated increases in storage requirements and computation time are expected. Stochastic Galerkin storage requirements scale linearly with the chosen basis order $P$ because the model stores $P+1$ expansion coefficients per variable per element. The ensemble averages $\left\langle\Phi_{p}\Phi_{s}\Phi_{l}\right\rangle$ in equation (23) and $\left\langle\Phi_{l}^{2}\right\rangle$ in equation (21) are constant values that can be precomputed once and stored. Due to the commutative property of the ensemble average of basis function products (equation 9), the associated storage requirements are small.

The expected increase in computation time can be estimated by examining the model formulation and assuming that calculations are performed without parallelisation. Choosing basis order $P$, the entire evolution equation (equation 21) must be evaluated $P+1$ times for the basis functions $0,\ldots,P$. For each evaluation, the ensemble average of numerical fluxes is calculated by sampling the Riemann solver $P+1$ times using Gauss-Hermite quadrature (equation 22). The computation time for the ensemble average of the source term vector (equation 23) can be neglected because its calculation is trivially fast compared to the Riemann solver. Hence, by considering the total samples of the Riemann solver, it is expected that the stochastic Galerkin model will be at least $\left(P+1\right)^{2}$ times slower than the deterministic solver.

Measuring the elapsed CPU time to simulate the steady-state test confirms this expectation: the stochastic Galerkin model with basis order $P=3$ is in fact about 20 times slower than the deterministic model, and the Riemann solver accounts for about 90% of the stochastic Galerkin total computation time. Compared to the Monte Carlo simulation which uses $2000$ deterministic iterations, the stochastic Galerkin model is about 100 times faster. Note that the deterministic model, stochastic Galerkin model and Monte Carlo simulation were implemented without parallelisation.

The speed-up observed here compares favourably with numerical tests performed by Ge et al. (2008). Their stochastic Galerkin model with basis order $P=5$ was about 200 times slower than their deterministic model. Their Monte Carlo simulation used $10\,000$ iterations, making their stochastic Galerkin model about 50 times faster. A direct comparison of computation time is not attempted since Ge et al. (2008) used a second-order finite volume formulation that is inherently more expensive than the first-order formulation presented here.

3.3.7 Opportunities for Parallel Computation

Computation time for Monte Carlo and stochastic Galerkin simulations could be reduced by exploiting opportunities for parallelism. Monte Carlo simulations are called ‘embarrassingly parallel’ because, given sufficient processors, it is easy to perform iterations entirely in parallel. The stochastic Galerkin model is not embarrassingly parallel, but most operations can be parallelised nevertheless. Due to basis orthogonality, the evolution equation (equation 21) can be evaluated in parallel over basis functions $\Phi_{0},\ldots,\Phi_{P}$. Gauss-Hermite quadrature of numerical fluxes (equation 22) can be parallelised, too. In this way, sampling the Riemann solver, which accounts for about 90% of the total computational cost, could be fully parallelised.

In theory, by fully exploiting parallelism, Monte Carlo and stochastic Galerkin simulations could be made to run almost as fast as a single iteration of the deterministic model, but this assumes access to sufficient processors. In practice, a fully parallelised Monte Carlo would require thousands of processors which are unavailable on typical hardware. In contrast, the stochastic Galerkin model with basis order $P=3$ could be fully parallelised with just 16 processors using current, commodity hardware. With more hardware, additional parallelism could be achieved using domain decomposition techniques. The stochastic Galerkin method imposes no barriers to domain decomposition because it preserves the local element-wise operations of the underlying deterministic formulation. Given these substantial reductions in computation time and hardware demands, a stochastic Galerkin shallow flow model could alleviate the computational constraints associated with conventional Monte Carlo simulations (Neal et al., 2013), and allow probabilistic simulations to account for more sources of uncertainty.

3.4 Test 3: Steady-State Slow Flow Over an Irregular Bed

While the previous tests are restricted to idealised topography profiles, this final test simulates slow flow over a highly irregular bed that is more representative of real-world river hydraulics. The bed is defined with a range of uncertainty that characterises bathymetric measurement errors. Stochastic Galerkin model results are verified against the analytic solution.

For this purpose, following Tseng (2004), the 1D domain is $[$0\text{,}\mathrm{m}$,$1500\text{,}\mathrm{m}$]$ with $M=200$ elements such that the mesh spacing is $\Delta x=$7.5\text{,}\mathrm{m}. The mean topography has an irregular profile as specified by Goutal and Maurel ([1997](#bib.bib8)) and is shown in figure [5](#S3.F5). The standard deviation of topography is $\sigma_{z}(x)=$0.5\text{\,}\mathrm{m} across the entire domain, which is chosen since it provides an acceptable margin of error for local-scale flood modelling applications (Schumann and Bates, 2018). The initial free-surface elevation is $15\text{\,}\mathrm{m}$ with no uncertainty. Subcritical boundary conditions are imposed such that the mean upstream discharge per unit-width is $0.75\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ and the mean downstream water depth is $15\text{\,}\mathrm{m}$, with no uncertainty on either boundary condition. Transmissive boundary conditions are used for the upstream water depth and downstream discharge. The test uses a timestep of $\Delta t=$0.5\text{,}\mathrm{s} and terminates at $t=$100\,000\text{\,}\mathrm{s} when the water has converged on a steady state down to a convergence error of ${10}^{-8}\text{\,}\mathrm{m}$ as measured by equation (44).

The stochastic Galerkin model is configured with a basis order $P=3$, and the resulting steady-state solution is shown in figure 5. Since the flow is very slow, the free-surface elevation is uniform. The free-surface elevation has the same standard deviation as the topography, which is due to the water depth at the downstream boundary being fixed at $15\text{\,}\mathrm{m}$ with no uncertainty. As a result, the water depth profile is certain, and the uncertainty in the free-surface elevation is determined entirely by the topography. The corresponding depth-averaged velocity profile is shown in figure 5. The stochastic Galerkin model calculates a velocity profile that is in excellent agreement with the analytic solution. The predicted velocity profile is certain with negligible standard deviation. This is expected since the water depth profile is unaffected by the prescribed topographic uncertainty. These results verify the robustness of the stochastic Galerkin model for flow over a highly-irregular, uncertain bed.

4 Conclusions and outlook

This paper proposed a stochastic Galerkin shallow flow model using the surface gradient method to ensure well-balancing. The well-balancedness of the stochastic model was studied theoretically and verified in a series of idealised numerical tests, with a prescribed uncertainty that characterises topographic or bathymetric measurement errors. These tests demonstrated that the stochastic Galerkin model preserves discrete balance between flux and topography gradients to machine precision, and accurately reproduces steady flow over highly irregular and uncertain topography. Despite its relatively low-order Wiener-Hermite basis, the stochastic Galerkin model was able to simulate strongly nonlinear behaviour, producing discontinuous and highly non-Gaussian probability distributions that compared favourably against a Monte Carlo reference solution.

It is likely that a more sophisticated treatment using a multiresolution wavelet decomposition (Le Maître et al., 2004; Pettersson et al., 2014) or multi-element discretisation of the stochastic dimension (Wan and Karniadakis, 2006; Li and Stinis, 2015) would further improve the representation of highly non-Gaussian and discontinuous probability distributions. It is also likely that a decomposition or discretisation of stochastic space would simplify a stochastic reformulation of positivity preserving schemes involving nonlinear operators, so enabling the development of stochastic wetting-and-drying processes that are essential for flood inundation modelling. Such improvements are the subject of future work.

The one-dimensional stochastic Galerkin model presented here is ideally suited to river routing applications: using as few as two expansion coefficients per element, the stochastic model can account for river bed measurement errors to produce a probabilistic flow simulation that is accurate and computationally efficient. Without parallelisation, the stochastic Galerkin model is at least 100 times faster than a Monte Carlo simulation. With parallelisation on commodity hardware, the computation time for the stochastic Galerkin model would be on par with a single iteration of the deterministic solver.

Acknowledgements

This work is part of the SEAMLESS-WAVE project (SoftwarE infrAstructure for Multi-purpose fLood modElling at variouS scaleS based on WAVElets) which is funded by the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/R007349/1. For information about the SEAMLESS-WAVE project visit https://www.seamlesswave.com. The authors are grateful to the editors and two anonymous reviewers for their helpful comments.

Appendix A Instructions for running the shallow flow models

Python 3 software is available for download from Zenodo (Shaw and Kesserwani, 2019b) which includes

a one-dimensional shallow flow model which can operate as a stochastic Galerkin model or a deterministic model 2. 2.

code for running Monte Carlo iterations of the deterministic model 3. 3.

preconfigured lake-at-rest, steady-state critical flow and steady-state flow over an irregular bed test cases

First ensure that Python 3, NumPy and SciPy are installed. Then, install the shallow flow model:

pip3 install --user --editable swepc.python

A.1 Running the Stochastic Galerkin Model

To run the stochastic Galerkin model:

swepc --degree <P> <testCase> <discretisation> --output-dir *<directory>*where <testCase> is either lakeAtRest, criticalSteadyState or tsengSteadyState, <discretisation> is either wellBalancedH for the well-balanced surface gradient method, or centredDifferenceH for the centred difference method given by equations (42) and (43). If --degree <P> is omitted then the default basis order $P=3$ is used. If --degree 0 is specified then the model operates as a deterministic model. On running the model, the following text files are written to the specified output directory:

coefficients.dat

Stochastic Galerkin expansion coefficients for topography $z$, water depth $h$ and discharge $q$

statistics.dat

Mean, standard deviation, skew and kurtosis calculated using equations (12), (13) and (14)

derived-statistics.dat

Statistics for the depth-averaged velocity and the free-surface elevation, a derived quantity calculated from equation (18)

Each file has one line per element with explanatory header rows prefixed by #.

A.2 Calculating Probability Density Functions

A probability density function for any variable can be calculated from the expansion coefficients for a given element:

sed -n *<line>*p coefficients.dat   | swepdf --min <value> --max <value> <variable> > pdf.datwhere <line> is the line number in coefficients.dat corresponding to the chosen element, and <variable> is z (topography), water (water depth), q (discharge) or derived-eta (free-surface elevation). The probability density function is calculated between the --min and --max values and written to pdf.dat.

A.3 Running Monte Carlo Simulations

To run a Monte Carlo simulation:

swepc --monte-carlo --mc-iterations <value>  criticalSteadyState <discretisation> --output-dir *<directory>*which writes the following text files:

statistics.dat

Mean, standard deviation, skew and kurtosis for topography, water depth and discharge

derived-statistics.dat

Mean and standard deviation for free-surface elevation and depth-averaged velocity \cprotect

**missing****ommandchars=

{}**

+sample*<n>*.dat+]Deterministic model data with one file per element (numbered $0,\ldots,M-1$), each having one line per Monte Carlo iteration