# Hyperbolicity-Preserving and Well-Balanced Stochastic Galerkin Method for
Two-Dimensional Shallow Water Equations

Dihan Dai11footnotemark: 1<EMAIL_ADDRESS>Yekaterina
Epshteyn22footnotemark: 2<EMAIL_ADDRESS>Akil Narayan33footnotemark:
344footnotemark: 4<EMAIL_ADDRESS>Department of Mathematics, University of
Utah, Salt Lake City, UT 84112, USA  Scientific Computing and Imaging (SCI)
Institute, University of Utah, Salt Lake City, UT 84112, USA

###### Abstract

Stochastic Galerkin formulations of the two-dimensional shallow water systems
parameterized with random variables may lose hyperbolicity, and hence change
the nature of the original model. In this work, we present a hyperbolicity-
preserving stochastic Galerkin formulation by carefully selecting the
polynomial chaos approximations to the nonlinear terms of
$(q^{x})^{2}/h,(q^{y})^{2}/h$, and $(q^{x}q^{y})/h$ in the shallow water
equations. We derive a sufficient condition to preserve the hyperbolicity of
the stochastic Galerkin system which requires only a finite collection of
positivity conditions on the stochastic water height at selected quadrature
points in parameter space. Based on our theoretical results for the stochastic
Galerkin formulation, we develop a corresponding well-balanced hyperbolicity-
preserving central-upwind scheme. We demonstrate the accuracy and the
robustness of the new scheme on several challenging numerical tests.

###### keywords:

finite volume method, stochastic Galerkin method, shallow water equations,
hyperbolic systems of conservation and balance laws
AMS subject classification. 35L65, 35Q35, 35R60, 65M60, 65M70

##  1 Introduction

The system of shallow water equations and closely related models have found
tremendous success in modelling flows whose vertical scales are significantly
smaller than their horizontal scales. Such flows include, but are not limited
to, the water flows in rivers, lakes, and coastal areas. In addition,
uncertainties can enter models, for instance, through measuring or empirically
approximating the bottom topography or initial conditions. This makes
deterministic-model-based predictions less reliable and therefore it is also
important to model this uncertainty in the governing shallow water systems. In
this work, we will consider the _parameterized_ stochastic shallow water
equations (SWE), where the parameter is a random variable that models the
uncertainty.

For discretizing parameterized partial differential equations, the polynomial
chaos expansion (PCE) technique is a commonly-used approach. It seeks to
represent a second-order random field using an expansion in an orthogonal
polynomial basis. It was first introduced by Wiener [44] to represent the
Gaussian processes using Hermite polynomials. The idea has been later extended
to the representations of random outputs and fields with more general
distributions using more general families of orthogonal polynomials [15, 49,
28, 43]. With PCE the strategy for encoding uncertainty, this gives rise to
parameterized partial differential equations. Two major classes of methods are
used to solve the resulting system and compute the statistics of the outputs.
One class of _non-intrusive_ methods is based on sampling (see, e.g.,[48, 33,
31]). The implementations and the parallelism can be easily achieved for non-
intrusive methods. However, such approximations can be less accurate than the
solutions of intrusive methods.

The class of _intrusive methods_ , in the PCE context, usually refers to
stochastic Galerkin (SG) method. In the SG framework, the truncated PCE
representations are projected by a Galerkin method in stochastic (parameter)
space to obtain a new system of deterministic partial differential equations
(PDEs) that describes the evolution of the coefficients in the PCE
representation. The new system typically requires new solvers and a
substantial rewrite of legacy codes. However, it is expected that the SG
method will deliver a more accurate solution than, for example, some
alternative non-intrusive strategies. In particular, it is known that
intrusive methods have near-optimal accuracy in an $L^{2}$ sense for static
problems [1, 29]. In addition, for SG formulations, one can generally prove
important properties of the system or of their subsequent numerical
discretizations. Because of this, applications of SG methods have found
success in modeling uncertainty in diffusion equations, e.g. [50, 11], kinetic
equations, e.g. [17, 40], and for conservation and balance laws with symmetric
Jacobian matrices, e.g. [42].

However, when the SG method is applied to a general hyperbolic system of
conservation/balance laws such as the SWE, the associated SG system may lose
hyperbolicity [10, 18, 14]. The loss of hyperbolicity may lead to unphysical
solutions and prevent the use of robust numerical schemes. Efforts have been
made to overcome this difficulty. In [45], a strategy to regain hyperbolicity
for the quasilinear form of hyperbolic systems is proposed. Hyperbolicity can
be regained by multiplying the SG formulation of the system by the left
eigenvector matrix of the flux Jacobian matrix, but this approach is limited
to quasilinear forms, and the numerical schemes designed for conservative
forms cannot be applied directly. Operator-splitting approaches for one-
dimensional Euler equations [7] and one-dimensional shallow water equations
[6] have also been investigated. The original systems are split into linear
hyperbolic and scalar subsystems, whose SG formulations remain hyperbolic.
However, it has been shown in [37] that the splitting approach may lead to an
overall non-hyperbolic system. Strategies that introduce a PCE of auxiliary
variables have also been explored. For example, the SG system of a
balance/conservation laws in terms of the entropic variables can be shown to
be hyperbolic [36, 35]. An optimization-based intrusive polynomial moment
method (IPMM) was proposed to compute the PCE of entropic variables given the
PCE of the conserved variables [10, 36, 35]. But the IPMM method can be
computationally expensive since one must solve an optimization problem for
each cell at each discrete time level. The PCE using a Roe variable
formulation is introduced in [34, 14, 13], where the flux of the SG system is
constructed using PCE of the Roe variables. Thus, the conservative form of the
system is preserved. It has been shown that, for the Wiener-Haar expansion,
both the SG Roe formulations of the Euler equation [34] and of the 1D SWE [14]
are hyperbolic. For the system of isothermal Euler equations, its SG Roe
formulation is hyperbolic under a positive definiteness condition of a matrix
depending on the positive square root of the water height, regardless of the
choices of distributions of the random parameters [14]. However, the Roe
formulation may still be expensive to implement since the PCE of the Roe
variables need to be calculated by solving both a nonlinear equation and a
linear equation. More discussion about the SG Roe formulation can be found in
[13]. Recently we developed a hyperbolicity-preserving SG formulation for the
one-dimensional shallow water equations [8], where we carefully tailored PCE
truncation operators for the nonlinear non-polynomial term in the system using
only the conserved variables. For arbitrary choices of distributions of the
random parameters, the resulting SG system is hyperbolic under only a
positive-definiteness condition on a matrix which depends on the PCE of the
water height. The condition is a stochastic variant of a positivity condition
for the water height for the deterministic system. The current paper is an
extension of the work [8] to the two dimensional shallow water equations.

In this work, the numerical scheme for the SG formulation in two-dimensional
physical space is based on the extension of the finite-volume central-upwind
scheme. The semi-discrete central-upwind schemes [22] are a class of robust
Godunov-type Riemann-solver-free projection-evolution methods for hyperbolic
systems. They are the generalizations of the central Nessyahu-Tadmor scheme
[32] and Kurganov-Tadmor high resolution central scheme [25]. The central
Nessyahu-Tadmor scheme, the Kurganov-Tadmor central scheme, and the family of
the central-upwind schemes have been successfully applied for the solution of
the general multi-dimensional hyperbolic systems of conservation/balance laws,
including for the solution of the deterministic SWE and related models [19]. A
second-order central-upwind scheme for the solution of SWE was first proposed
in [20], but it was not designed to preserve the positivity and the well-
balanced property simultaneously. The scheme was further improved in [23] to
capture the “lake-at-rest” steady state and to preserve the positivity of the
water height. See [24, 21, 4, 5, 30, 19] for more examples and details of
other closely related work. The scheme in this paper is mainly based on
further extensions to stochastic SWE of the schemes proposed in [22, 23]. We
note that, although the developed framework in [8] and in the current paper is
not limited to the choice of central-upwind schemes, the Riemann-solver-free
nature of these schemes allows us to design an efficient deterministic solver
as a part of the numerical discretization of the random SWEs.

###  1.1 Contributions of This Paper

In this paper we extend the ideas that we developed in [8] to stochastic
shallow water equations in two-dimensional physical space. In order to develop
a hyperbolicity-preserving stochastic Galerkin formulation, we use specially
tailored PCE truncation operators for the SWE nonlinear term $(q^{x}q^{y})/h$
that occurs in both the $x$\- and the $y$\- component of the flux. For any
probability distributions of the random parameters, our SG formulation is
conservative and we prove it is also hyperbolic under a positive-definiteness
condition on the matrix associated with the PCE of the water height (this
result is similar to the one-dimensional physical space result in [8]). This
positive-definiteness condition is a stochastic variant of the positivity
condition on the water height for the deterministic model. We show that the
positive-definiteness condition can be guaranteed by enforcing a finite
collection of positivity conditions on the stochastic water height at the
stochastic quadrature nodes. Based on our SG formulation, we also develop a
corresponding well-balanced hyperbolicity-preserving central-upwind scheme.

This paper is structured as follows. In Section 2, we introduce the
parameterized stochastic shallow water model in 2D physical space, as well as
the stochastic Galerkin discretization of the system using a particular choice
of the PCE truncation for the nonlinear terms $(q^{x})^{2}/h,(q^{y})^{2}/h$,
and $(q^{x}q^{y})/h$. In Section 3, we discuss the hyperbolicity of the
resulting stochastic Galerkin system and present a sufficient condition to
guarantee hyperbolicity of the stochastic Galerkin formulation of the 2D
shallow water system. In Section 4 we develop a well-balanced central-upwind
scheme for the resulting model and obtain a hyperbolicity-preserving CFL-type
condition. In Section 5 we demonstrate the accuracy and the robustness of the
developed numerical scheme on several challenging numerical tests. Conclusions
and future directions are discussed in Section 6.

##  2 Modeling Stochastic Shallow Water Equations

###  2.1 Stochastic Modelling of the Two-Dimensional Shallow Water Equations

The two-dimensional system of shallow water equations read

$\displaystyle\dfrac{\partial h}{\partial t}+\dfrac{\partial q^{x}}{\partial
x}+\dfrac{\partial q^{y}}{\partial y}=0,$ (1) $\displaystyle\dfrac{\partial
q^{x}}{\partial t}+\dfrac{\partial}{\partial
x}\left(\frac{(q^{x})^{2}}{h}+\frac{gh^{2}}{2}\right)+\dfrac{\partial}{\partial
y}\left(\frac{q^{x}q^{y}}{h}\right)=-gh\dfrac{\partial B}{\partial x},$
$\displaystyle\dfrac{\partial q^{y}}{\partial t}+\dfrac{\partial}{\partial
x}\left(\frac{q^{x}q^{y}}{h}\right)+\dfrac{\partial}{\partial
y}\left(\frac{(q^{y})^{2}}{h}+\frac{gh^{2}}{2}\right)=-gh\dfrac{\partial
B}{\partial y}.$

where $h$ is the height of the water, $g$ is the gravitational constant, $B$
is the time-independent bottom topography function, and $q^{x}$ and $q^{x}$
are the discharges in the $x$\- and $y$-directions, respectively. The
quantities $h$, $q^{x}$, and $q^{y}$ depend on the spatial variables $(x,y)$
and time $t$, and $B$ depends only on $(x,y)$. We are interested in modeling a
random/uncertain variant of this system.

Let $(\Omega,\mathcal{F},P)$ be a complete probability space with event space
$\Omega$, $\sigma$-algebra $\mathcal{F}$, and probability measure $P$ defined
on $\mathcal{F}$. Uncertainty can enter the deterministic system, for example,
via the bottom topography $B$. For example, the stochastic bottom can be
modeled by a $d$-dimensional random field

$\displaystyle
B=B(x,y,\xi(\omega))=B_{0}(x,y)+\sum_{k=1}^{d}B_{k}(x,y)\xi_{k},$

where $\xi(\omega)=(\xi_{1}(\omega),\ldots,\xi_{d}(\omega))$ is a
$d$-dimensional random variable for $\omega\in\Omega$. Such a stochastic
bottom function will lead to a $\xi$-parameterized version of the 2D shallow
water equations (1), so that now $h=h(x,y,t,\xi)$, $q^{x}=q^{x}(x,y,t,\xi)$,
$q^{y}=q^{x}(x,y,t,\xi)$, and $B=B(x,y,\xi)$. Our results can be generalized
to other models of uncertainty, e.g., the uncertainty from the measurements of
the initial water displacements $\eta\coloneqq h+B$. In the following context,
we will refer to the quantities associated only with $(x,y)$ as “physical”
quantities and the quantities associated only with $\xi$ as “stochastic”
quantities.

####  2.1.1 Polynomial Chaos Expansion - A Review

In this section, we briefly review the results and notation for polynomial
chaos expansions (PCE), which are similar to what is used in [8]. We include
this section here in order to be self-contained; interested readers can find a
more technical discussin of PCE in [9, 47, 41].

Let $\rho$ be the Lebesgue density function for the random variable
$\xi\in\mathbbm{R}^{d}$. The strategy of PCE is to approximate the dependence
of a square-integrable random field by a polynomial of $\xi$. Let
$\nu=(\nu_{1},\ldots,\nu_{d})\in\mathbbm{N}_{0}^{d}$ be a multi-index. For a
$d$-dimensional vector $\zeta\in\mathbbm{R}^{d}$, the standard notation for
the monomials with multi-index $\nu$ reads

$\displaystyle\zeta^{\nu}$
$\displaystyle\coloneqq\prod_{j=1}^{d}\zeta_{j}^{\nu_{j}},$
$\displaystyle\zeta^{0}=\zeta^{(0,0,\ldots,0)}$ $\displaystyle=1.$

We consider the PCE approximations in the polynomial subspace defined by
$\Lambda$:

$\displaystyle P_{\Lambda}$
$\displaystyle=\mathrm{span}\\{\zeta^{\nu}\;\;\big{|}\;\;\nu\in\Lambda\\},$
$\displaystyle\dim P_{\Lambda}$ $\displaystyle=K\coloneqq|\Lambda|,$

where $\Lambda\subset\mathbbm{N}_{0}^{d}$ denotes any non-empty, size-$K$
finite set of multi-indices with $0=(0,0,\ldots,0)\in\Lambda$. We will also
need “powers” of this set, defined by $r$-fold products of $P_{\Lambda}$
elements:

$\displaystyle P_{\Lambda}^{r}$
$\displaystyle\coloneqq\mathrm{span}\left\\{\prod_{j=1}^{r}p_{j}\;\;\big{|}\;\;p_{j}\in
P_{\Lambda},\;j=1,\ldots,r\right\\},$ $\displaystyle\dim P_{\Lambda}^{r}$
$\displaystyle\leq\left(\\!\\!\left(\begin{array}[]{c}K\\\
r\end{array}\right)\\!\\!\right)=\left(\begin{array}[]{c}K+r-1\\\
r\end{array}\right),$ (6)

where the dimension bound results from a combinatoric argument. Note that
since $0\in\Lambda$, then $P^{r}_{\Lambda}\subseteq P^{s}_{\Lambda}$ for any
$r\leq s$. We will later exercise the notation above for $r=3$ (3.2). Define
the function space

$\displaystyle
L^{2}_{\rho}(\mathbbm{R}^{d})\coloneqq\left\\{f:\mathbbm{R}^{d}\to\mathbbm{R}\middle|\left(\int_{\mathbbm{R}^{d}}f^{2}(s)\rho(s)ds\right)^{\frac{1}{2}}<+\infty\right\\}.$

Assuming finite polynomial moments of all orders for $\rho$, there exists an
$L^{2}_{\rho}(\mathbbm{R}^{d})$-orthonormal polynomial basis
$\\{\phi_{k}\\}_{k=1}^{\infty}$ of $P_{\Lambda}$, i.e.,

$\displaystyle\langle\phi_{k},\phi_{\ell}\rangle_{\rho}$
$\displaystyle\coloneqq\int_{\mathbbm{R}^{d}}\phi_{k}(s)\phi_{\ell}(s){\rho}(s)ds=\delta_{k,\ell},$
$\displaystyle\phi_{1}(\xi)$ $\displaystyle\equiv 1,$ (7)

for all $k,\ell\in\mathbbm{N}$, where $\delta_{k,\ell}$ is the Kronecker
delta.

If $z(x,y,t,\cdot)\in L^{2}_{\rho}(\mathbbm{R}^{d})$, then under mild
conditions on the probability measure $\rho$ (see [12]) there exists a
convergent expansion of $z$ in these basis functions,

$\displaystyle z(x,y,t,\xi)$ $\displaystyle\stackrel{{\scriptstyle
L^{2}_{\rho}}}{{=}}\sum_{k=1}^{\infty}\widehat{z}_{k}(x,y,t)\phi_{k}(\xi),$

where $\widehat{z}_{k}(x,y,t)$ are (deterministic) Fourier coefficients in the
basis $\\{\phi_{k}\\}_{k\in\mathbbm{N}}$. If
$P_{\Lambda}=\mathrm{span}\\{\phi_{k}\;\;\big{|}\;\;1\leq k\leq K\\}$, then we
define the $K$-term PCE approximation of $z$ on $P_{\Lambda}$ to be the
truncation of the first $K$ terms in the summation above,

$\displaystyle
z(x,y,t,\xi)\approx\sum_{k=1}^{K}\widehat{z}_{k}(x,y,t)\phi_{k}(\xi)\eqqcolon\Pi_{\Lambda}[z](x,y,t,\xi).$
(8)

The operator $\Pi_{\Lambda}:L^{2}_{\rho}(\mathbbm{R}^{d})\to P_{\Lambda}$ is
the linear projection operator onto $P_{\Lambda}$.

Orthogonality of the basis functions allows us to compute statistics of the
the truncated PCE by manipulating the coefficients in a straightforward way.
For example, the expectation and the variance of $\Pi_{\Lambda}[z](x,y,t,\xi)$
are given by

$\mathbb{E}[\Pi_{\Lambda}[z](x,y,t,\xi)]=\widehat{z}_{1}(x,y,t),\quad\text{Var}[\Pi_{\Lambda}[z](x,y,t,\xi)]=\sum_{k=2}^{K}\widehat{z}_{k}^{2}(x,y,t),$
(9)

respectively.

We introduce the _pseudo-spectral product_ as the approximation to the
$K$-term PCE of the product of two random fields $a(x,y,t,\xi)$ and
$b(x,y,t,\xi)$,

$\displaystyle\Pi_{\Lambda}[a,b]$
$\displaystyle\coloneqq\Pi_{\Lambda}\left[\Pi_{\Lambda}[a]\;\Pi_{\Lambda}[b]\right]=\sum_{m=1}^{K}\left(\sum_{k,\ell=1}^{K}\widehat{a}_{k}\widehat{b}_{\ell}\langle\phi_{k}\phi_{\ell},\phi_{m}\rangle_{\rho}\right)\phi_{m}(\xi).$
(10)

It is an exact projection onto $P_{\Lambda}$ of the product of two
$P_{\Lambda}$ projections. The approximation (10) is a widely used strategy
for computing PCE expansion of products (e.g. [9, 14]). Such an operation can
be cast in linear algebraic terms by considering vectors comprised of the PCE
expansion coefficients. Let $\widehat{z}\in\mathbbm{R}^{K}$ denote the vector
of the $\phi_{k}$-expansion coefficients of $z\in P_{\Lambda}$. Define the
linear operator $\mathcal{P}:\mathbbm{R}^{K}\rightarrow\mathbbm{R}^{K\times
K}$ as

$\displaystyle\mathcal{P}(\widehat{z})$
$\displaystyle\coloneqq\sum_{k=1}^{K}\widehat{z}_{k}\mathcal{M}_{k},$
$\displaystyle\mathcal{M}_{k}$ $\displaystyle\in\mathbbm{R}^{K\times K},$
$\displaystyle(\mathcal{M}_{k})_{\ell m}$
$\displaystyle=\langle\phi_{k},\phi_{\ell}\phi_{m}\rangle_{\rho},$ (11)

where $\mathcal{M}_{k}$ is a symmetric matrix for each $k$. The following
properties can be verified:

$\displaystyle\mathcal{P}(\widehat{a})$
$\displaystyle=\begin{pmatrix}\mathcal{M}_{1}\widehat{a}|\mathcal{M}_{2}\widehat{a}|\ldots|\mathcal{M}_{K}\widehat{a}\end{pmatrix},$
$\displaystyle\mathcal{P}(\widehat{a})\widehat{b}$
$\displaystyle=\mathcal{P}(\widehat{b})\widehat{a},$
$\displaystyle\widehat{\Pi_{\Lambda}[a,b]}=\mathcal{P}(\widehat{a})\widehat{b},$
(12)

where the last property is due to (10).

We will also need to compute $P_{\Lambda}$ truncations of ratios of random
fields (when for each $(x,y,t)$ the denominator is a single-signed variable
with probability 1). For a single-signed random field $a$, we have

$\Pi_{\Lambda}\left[a\,\frac{b}{a}\right](x,y,t,\xi)=\Pi_{\Lambda}[b](x,y,t,\xi).$
(13)

Using (13) as motivation, we make the assumption,

$\displaystyle\Pi_{\Lambda}\left[a,\frac{b}{a}\right]=\Pi_{\Lambda}[b]\hskip
5.0pt\stackrel{{\scriptstyle\eqref{eq:pmatrixproperty}}}{{\Longleftrightarrow}}\hskip
5.0pt\mathcal{P}(\widehat{a})\widehat{\left(\frac{b}{a}\right)}=\widehat{b},$
(14)

which then motivates the definition for a new operator for the $K$-term PCE of
the ratios,

$\Pi^{\dagger}_{\Lambda}\left[\frac{b}{a}\right](\xi)\coloneqq\sum_{k=1}^{K}c_{k}\phi_{k}(\xi),$
(15)

where $c_{k}$ is the $k$th element of $\widehat{\left(\frac{b}{a}\right)}$
defined by (14), assuming $\mathcal{P}(\widehat{a})$ is invertible.

###  2.2 Stochastic Galerkin Formulation for 2D SWE

Following a standard Galerkin procedure in stochastic ($\xi$) space, we first
reduce the problem to a finite-dimensional one by making the ansatz that the
solutions to $h,q^{x},$ and $q^{y}$ lie in the polynomial space $P_{\Lambda}$,

$\displaystyle h\simeq
h_{\Lambda}\;\coloneqq\sum_{k=1}^{K}\widehat{h}_{k}(x,y,t)\phi_{k}(\xi),$ (16)
$\displaystyle q^{x}\simeq
q^{x}_{\Lambda}\coloneqq\sum_{k=1}^{K}(\widehat{q^{x}})_{k}(x,y,t)\phi_{k}(\xi),$
$\displaystyle q^{y}\simeq
q^{y}_{\Lambda}\coloneqq\sum_{k=1}^{K}(\widehat{q^{y}})_{k}(x,y,t)\phi_{k}(\xi),$

respectively, and $B$ is replaced by $\Pi_{\Lambda}[B]$. We then apply the
projection operator $\Pi_{\Lambda}$ to both sides of (1) and enforce equality.
Similar to [8], we make the assumption that,

$\displaystyle\frac{(q^{x})^{2}}{h}=\frac{q^{x}}{h}\;q^{x}\hskip
10.0pt\longrightarrow\hskip
10.0pt\Pi_{\Lambda}\left[\frac{(q^{x}_{\Lambda})^{2}}{h_{\Lambda}}\right]=\Pi_{\Lambda}\left[q^{x}_{\Lambda}\;\Pi^{\dagger}_{\Lambda}\left[\frac{q^{x}_{\Lambda}}{h_{\Lambda}}\right]\right],$
(17) $\displaystyle\frac{(q^{y})^{2}}{h}=\frac{q^{y}}{h}\;q^{y}\hskip
10.0pt\longrightarrow\hskip
10.0pt\Pi_{\Lambda}\left[\frac{(q^{y}_{\Lambda})^{2}}{h_{\Lambda}}\right]=\Pi_{\Lambda}\left[q^{y}_{\Lambda}\;\Pi^{\dagger}_{\Lambda}\left[\frac{q^{y}_{\Lambda}}{h_{\Lambda}}\right]\right].$

In addition, we make another crucial assumption on the approximation of the
nonlinear term $(q^{x}q^{y})/h$ that occurs in the fluxes of both $x$\- and
$y$-directions. Specifically, for the nonlinear term $q^{x}q^{y}/h$ in
$(q^{x}q^{y}/h)_{x}$, we use the approximation,

$\frac{q^{x}q^{y}}{h}=(q^{x})\frac{q^{y}}{h}\;\hskip
10.0pt\longrightarrow\hskip
10.0pt\Pi_{\Lambda}\left[\frac{q^{x}_{\Lambda}q^{y}_{\Lambda}}{h_{\Lambda}}\right]=\Pi_{\Lambda}\left[q^{x}_{\Lambda}\;\Pi^{\dagger}_{\Lambda}\left[\frac{q^{y}_{\Lambda}}{h_{\Lambda}}\right]\right],$
(18)

while for the same nonlinear term in $(q^{x}q^{y}/h)_{y}$, we adopt the
approximation,

$\frac{q^{x}q^{y}}{h}=(q^{y})\frac{q^{x}}{h}\;\hskip
10.0pt\longrightarrow\hskip
10.0pt\Pi_{\Lambda}\left[\frac{q^{x}_{\Lambda}q^{y}_{\Lambda}}{h_{\Lambda}}\right]=\Pi_{\Lambda}\left[q^{y}_{\Lambda}\;\Pi^{\dagger}_{\Lambda}\left[\frac{q^{x}_{\Lambda}}{h_{\Lambda}}\right]\right].$
(19)

The two different approximations (18)-(19) turn out to be the crucial
ingredient of our hyperbolicity-preserving SG formulation (3.1). The two
approximations are in general not equal. But as $K$ increases, we observe
empirically that the discrepancy decreases exponentially (see Figure 1).

Figure 1: Numerical verification: the two approximations (18)-(19) are in
general not the same. The difference is computed as the maximum of
$L^{2}_{\rho}$ norm of the difference between the two PCE approximations of
$q^{x}q^{y}/h$ in corresponding cells over the whole computational domain. The
data are taken from the numerical solutions to (62)-(63) at $t=0.07$ for
different 2D grid sizes.

The steps introduced above lead to a new system of balance laws for the PCE
vector of the stochastic state variables,

$\dfrac{\partial}{\partial t}(\widehat{U})+\dfrac{\partial}{\partial
x}(\widehat{F}(\widehat{U}))+\dfrac{\partial}{\partial
y}(\widehat{G}(\widehat{U}))=\widehat{S}(\widehat{U},\widehat{B}).$ (20)

Here
$\widehat{U}\coloneqq(\widehat{h}^{\top},\widehat{q^{x}}^{\top},\widehat{q^{y}}^{\top})^{\top}$,
where $\widehat{h},\widehat{q^{x}}$, and $\widehat{q^{y}}$ are each length-$K$
vectors whose entries are the coefficients introduced in (16). The flux terms
are,

$\displaystyle\widehat{F}(\widehat{U})=\begin{pmatrix}\widehat{q^{x}}\\\
\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\widehat{q^{x}}+\frac{1}{2}g\mathcal{P}(\widehat{h})\widehat{h}\\\
\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\widehat{q^{y}}\end{pmatrix},$
$\displaystyle\widehat{G}(\widehat{U})=\begin{pmatrix}\widehat{q^{y}}\\\
\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})\widehat{q^{x}}\\\
\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})\widehat{q^{y}}+\frac{1}{2}g\mathcal{P}(\widehat{h})\widehat{h}\end{pmatrix},$
(21)

and the source term is

$\displaystyle\widehat{S}(\widehat{U},\widehat{B})=\begin{pmatrix}0\\\
-g\mathcal{P}(\widehat{h})\widehat{B_{x}}\\\
-g\mathcal{P}(\widehat{h})\widehat{B_{y}}\end{pmatrix}.$ (22)

Following a similar calculation as in [18], we obtain the flux Jacobians
written in terms of $K\times K$ blocks,

$\displaystyle\frac{\partial\widehat{F}}{\partial\widehat{U}}$
$\displaystyle=\begin{pmatrix}O&I&O\\\
g\mathcal{P}(\widehat{h})-\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\mathcal{P}(\widehat{u})&\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})+\mathcal{P}(\widehat{u})&O\\\
-\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\mathcal{P}(\widehat{v})&\mathcal{P}(\widehat{v})&\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\end{pmatrix},$
(23) $\displaystyle\frac{\partial\widehat{G}}{\partial\widehat{U}}$
$\displaystyle=\begin{pmatrix}O&O&I\\\
-\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})\mathcal{P}(\widehat{u})&\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})&\mathcal{P}(\widehat{u})\\\
g\mathcal{P}(\widehat{h})-\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})\mathcal{P}(\widehat{v})&O&\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})+\mathcal{P}(\widehat{v})\end{pmatrix},$
(24)

where we have introduced

$\displaystyle\widehat{u}=\mathcal{P}^{-1}(\widehat{h})\widehat{q^{x}},$
$\displaystyle\widehat{v}=\mathcal{P}^{-1}(\widehat{h})\widehat{q^{y}}$ (25)

which can be interpreted as PCE coefficients of the $x$\- and $y$-velocity,
respectively.

###### Remark 2.1

In the numerical section Section 4, our scheme requires us to compute the
maximum and minimum eigenvalues of both
$\frac{\partial\widehat{F}}{\partial\widehat{U}}$ and
$\frac{\partial\widehat{G}}{\partial\widehat{U}}$. Some algebraic
manipulations can accelerate the computation. For example, the flux Jacobian
$\frac{\partial\widehat{F}}{\partial\widehat{U}}$ has the same eigenvalues as
the block diagonal matrix

$\begin{pmatrix}O&I&O\\\
g\mathcal{P}(\widehat{h})-\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\mathcal{P}(\widehat{u})&\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})+\mathcal{P}(\widehat{u})&O\\\
O&O&\mathcal{P}(\widehat{q^{x}})\mathcal{P}^{-1}(\widehat{h})\end{pmatrix}.$

One can first compute the maximum and minimum eigenvalues of the first
$2K\times 2K$ submatrix and the $K\times K$ submatrix on the right corner, and
then determine the maximum and minimum eigenvalues of the original matrix by
comparing the eigenvalues of the submatrices. Similar calculations can be
applied to $\frac{\partial\widehat{G}}{\partial\widehat{U}}$, where we can
replace the
$-\mathcal{P}(\widehat{q^{y}})\mathcal{P}^{-1}(\widehat{h})\mathcal{P}(\widehat{u})$
and $\mathcal{P}(\widehat{u})$ blocks by zero matrices when computing
eigenvalues.

##  3 Hyperbolicity of the SG System

In this section we will present a sufficient condition to ensure hyperbolicity
of the system (20) as well as a computationally tractable way to preserve
hyperbolicity of the system.

###### Theorem 3.1

If the matrix $\mathcal{P}(\widehat{h})$ is strictly positive definite for
every point $(x,y,t)$ in the computational spatial-temporal domain, then the
SG formulation (20) is hyperbolic.

###### Proof 3.1

To show that the system (20) is hyperbolic, we need to show that, for any unit
vector $\boldsymbol{n}=(n_{x},n_{y})^{\top}$, the matrix

$\displaystyle
n_{x}\frac{\partial\widehat{F}}{\partial\widehat{U}}+n_{y}\frac{\partial\widehat{G}}{\partial\widehat{U}}$
(26)

is real diagonalizable. Define,

$\displaystyle A\coloneqq g\mathcal{P}(\widehat{q^{x}}),$ $\displaystyle
B\coloneqq\mathcal{P}(\widehat{u}),$ $\displaystyle C\coloneqq
g\mathcal{P}(\widehat{q^{y}}),$ $\displaystyle
D\coloneqq\mathcal{P}(\widehat{v}),$ $\displaystyle
E\coloneqq\sqrt{g\mathcal{P}(\widehat{h})},$ (27)

where $\sqrt{M}$ is the (unique) symmetric positive definite square root of a
symmetric positive definite matrix $M$. Define

$P=\begin{pmatrix}I&O&I\\\ B+n_{x}E&-n_{y}E&B-n_{x}E\\\
D+n_{y}E&n_{x}E&D-n_{y}E\end{pmatrix},$

whose inverse is given by

$P^{-1}=\left(\frac{1}{2}\right)(I\otimes
E^{-1})\begin{pmatrix}E-n_{x}B-n_{y}D&n_{x}I&n_{y}I\\\
2(n_{y}B-n_{x}D)&-2n_{y}I&2n_{x}I\\\ E+n_{x}B+n_{y}D&-n_{x}I&-n_{y}I\\\
\end{pmatrix}.$

Using (23)-(24), direct calculation leads to

$J\coloneqq\begin{aligned}
P^{-1}\left(n_{x}\frac{\partial\widehat{F}}{\partial\widehat{U}}+n_{y}\frac{\partial\widehat{G}}{\partial\widehat{U}}\right)P=&\left(\frac{1}{2}\right)\begin{pmatrix}J_{11}&O&J_{13}\\\
O&J_{22}&O\\\ J_{31}&O&J_{33},\end{pmatrix}.\end{aligned}$ (28)

where,

$\displaystyle J_{11}=2E+n_{x}(B+E^{-1}AE^{-1})+n_{y}(D+E^{-1}CE^{-1}),$ (29)
$\displaystyle J_{13}=n_{x}(B-E^{-1}AE^{-1})+n_{y}(D-E^{-1}CE^{-1})=J_{31},$
$\displaystyle J_{22}=2n_{x}E^{-1}AE^{-1}+2n_{y}E^{-1}CE^{-1},$ $\displaystyle
J_{33}=-2E+n_{x}(B+E^{-1}AE^{-1})+n_{y}(D+E^{-1}CE^{-1}).$

From (27), the matrices $A,B,C,D$ and $E$ are all symmetric. Thus, it can be
verified that $J_{11},J_{13},J_{22},J_{31},$ and $J_{33}$ are all symmetric.
Since $J_{13}=J_{31}$, the matrix $J$ is symmetric and therefore the similar
matrix
$\left(n_{x}\frac{\partial\widehat{F}}{\partial\widehat{U}}+n_{y}\frac{\partial\widehat{G}}{\partial\widehat{U}}\right)$
is real diagonalizable.

###### Remark 3.1

When the system is deterministic, the matrix $P$ and the matrix $J$ in the
proof above reduce to the eigenmatrix that symmetrizes the deterministic
Jacobian matrix and the diagonal matrix after symmetrization, respectively.

From 3.1, the hyperbolicity of the SG system relies only on the positive
definiteness of the matrix $\mathcal{P}(\widehat{h})$, a stochastic variant of
the hyperbolic condition for the deterministic 2D SWE. In [8], a
computationally convenient sufficient condition to guarantee
$\mathcal{P}(\widehat{h})>0$ is given by the following theorem.

###### Theorem 3.2

Given $\Lambda$, let nodes $\xi_{m}$ and weights $\tau_{m}$ satisfying
$\\{(\xi_{m},\tau_{m})\\}_{m=1}^{M}\subset\mathbbm{R}^{d}\times(0,\infty)$
represent any $M$-point positive quadrature rule that is exact on
$P_{\Lambda}^{3}$, i.e.,

$\displaystyle\int_{\mathbbm{R}^{d}}p(\xi)\rho(\xi)d\xi$
$\displaystyle=\sum_{m=1}^{M}p(\xi_{m})\tau_{m},$ $\displaystyle p$
$\displaystyle\in P_{\Lambda}^{3}.$ (30)

If

$\displaystyle h_{\Lambda}(x,y,t,\xi_{m})>0\;\;\forall\;m=1,\ldots,M,$ (31)

then the SG system of the 2D SWE (20) is hyperbolic.

###### Proof 3.2

See Theorem 3.4 in [8].

Thus, to preserve the hyperbolicity of (20), we only need to preserve a finite
number of positivity conditions on the stochastic water height. In the worst-
case scenario, the number $M$ of quadrature points used in the stochastic
space is of order $O(\text{dim}P_{\Lambda}^{3})\lesssim K^{3}$ (see (6) for
the definition of $P_{\Lambda}^{3}$). For one-dimensional stochastic space,
one can use the $M$-point Gaussian quadrature associated with the density
$\rho$ of $\xi$, with $M$ any integer satisfying
$M\geq\left\lceil\frac{3K}{2}\right\rceil-1$.

Instead of preserving the positive definiteness of $\mathcal{P}(\widehat{h})$,
one could choose to preserve the positivity of the expectation
$\widehat{h}_{1}$ of $h_{\Lambda}$. The following corollary states that this
condition is guaranteed by $\mathcal{P}(\widehat{h})>0$.

###### Corollary 3.2.1

If the assumptions of Theorem (3.2) hold, then $\widehat{h}_{1}>0$.

###### Proof 3.3

See Corollary 3.6 in [8].

##  4 Numerical Scheme for Stochastic Shallow Water Equations

In this section, we present a well-balanced central-upwind scheme that
preserves the hyperbolicity of the system (20) at every time step.

###  4.1 Central-Upwind Scheme for the SG System

For simplicity, we consider only rectangular spatial domains. Let
$\left\\{\mathcal{C}_{i,j}\right\\}_{i\in[N_{x}],j\in[N_{y}]}$ be a uniform
rectangular partition over the rectangular computational domain, where

$\displaystyle\mathcal{C}_{i,j}\coloneqq\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]\times\left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right],$

$[N_{x}]$ is the index set $\\{1,2,\ldots,N_{x}\\}$, and $[N_{y}]$ is defined
similarly. Define $\Delta x\coloneqq
x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},\Delta
y=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$. Let $|\mathcal{C}_{i,j}|$ denote the
size $\Delta x\Delta y$ of cell $(i,j)$. Define the approximation to the cell
average of the vector $\widehat{U}$ over cell $\mathcal{C}_{i,j}$ to be

$\displaystyle\overline{\boldsymbol{U}}_{i,j}(t)\coloneqq\begin{pmatrix}\overline{\boldsymbol{h}}_{i,j}(t)\\\
\overline{\boldsymbol{q}^{x}}_{i,j}(t)\\\
\overline{\boldsymbol{q}^{y}}_{i,j}(t)\end{pmatrix}\approx\frac{1}{|\mathcal{C}_{i,j}|}\bigintsss_{\mathcal{C}_{i,j}}\begin{pmatrix}\widehat{h}(x,y,t)\\\
\widehat{q^{x}}(x,y,t)\\\
\widehat{q^{y}}(x,y,t)\end{pmatrix}dx\;dy\in\mathbbm{R}^{3K}.$

Since $\widehat{U}$ with subscripts has been used for the polynomial moments,
we use the bold letter $\boldsymbol{U}$ with superscripts and subscripts here
in the following subsections to represent the cell averages and the pointwise
values of the reconstructions, respectively, for the vector function
$\widehat{U}(x,y,t)$. For example, $\boldsymbol{U}^{W}_{i,j}$ is the
reconstructed value of $\widehat{U}$ at the west-side
$(x_{i-\frac{1}{2}},y_{j})$ of the cell $\mathcal{C}_{i,j}$. A similar
reasoning applies to the pairs $(\boldsymbol{h},\widehat{h})$,
$(\boldsymbol{q}^{x},\widehat{q^{x}})$, and
$(\boldsymbol{q}^{y},\widehat{q^{y}})$. To minimize clutter, we will
notationally suppress $t$ dependence from here onward.

For a continuous bottom topography, the pointwise values of the
reconstructions are given by

$\displaystyle\boldsymbol{U}^{E}_{i,j}=\overline{\boldsymbol{U}}_{i,j}+\frac{\Delta
x}{2}(\boldsymbol{U}_{x})_{i,j},$
$\displaystyle\boldsymbol{U}^{W}_{i,j}=\overline{\boldsymbol{U}}_{i,j}-\frac{\Delta
x}{2}(\boldsymbol{U}_{x})_{i,j},$
$\displaystyle\boldsymbol{U}^{N}_{i,j}=\overline{\boldsymbol{U}}_{i,j}+\frac{\Delta
y}{2}(\boldsymbol{U}_{y})_{i,j},$
$\displaystyle\boldsymbol{U}^{S}_{i,j}=\overline{\boldsymbol{U}}_{i,j}-\frac{\Delta
y}{2}(\boldsymbol{U}_{y})_{i,j},$

where $\boldsymbol{U}_{x}$ and $\boldsymbol{U}_{y}$ are computed using the
generalized minmod limiter,

$\displaystyle(\boldsymbol{U}_{x})_{i,j}=\text{minmod}\left(\theta\frac{\overline{\boldsymbol{U}}_{i,j}-\overline{\boldsymbol{U}}_{i-1,j}}{\Delta
x},\frac{\overline{\boldsymbol{U}}_{i+1,j}-\overline{\boldsymbol{U}}_{i-1,j}}{2\Delta
x},\theta\frac{\overline{\boldsymbol{U}}_{i+1,j}-\overline{\boldsymbol{U}}_{i,j}}{\Delta
x}\right),$
$\displaystyle(\boldsymbol{U}_{y})_{i,j}=\text{minmod}\left(\theta\frac{\overline{\boldsymbol{U}}_{i,j}-\overline{\boldsymbol{U}}_{i,j-1}}{\Delta
y},\frac{\overline{\boldsymbol{U}}_{i,j+1}-\overline{\boldsymbol{U}}_{i,j-1}}{2\Delta
y},\theta\frac{\overline{\boldsymbol{U}}_{i,j+1}-\overline{\boldsymbol{U}}_{i,j}}{\Delta
y}\right).$

The minmod function for scalars is defined by

$\text{minmod}(c_{1},c_{2},\ldots,c_{n})\coloneqq\left\\{\begin{aligned}
&\max\nolimits_{j}\\{c_{j}\\},&&\text{if }c_{j}<0\;\;\forall j,\\\
&\min\nolimits_{j}\\{c_{j}\\},&&\text{if }c_{j}>0\;\;\forall j,\\\
&0,&&\text{otherwise},\\\ \end{aligned}\right.$

and the minmod function for vector is simply applying the scalar minmod
function for each component. The parameter $\theta\in[1,2]$ controls how
dissipative the limiter is. The limiter is less dissipative for larger value
of $\theta$, but it remains non-oscillatory for any choice of $\theta$ [22].

The speeds of local propagation are then estimated by,

$\displaystyle
a^{-}_{i+\frac{1}{2},j}=\min\left\\{\lambda_{1}\left(\frac{\partial\widehat{F}}{\partial\widehat{U}}(\boldsymbol{U}^{W}_{i+1,j})\right),\lambda_{1}\left(\frac{\partial\widehat{F}}{\partial\widehat{U}}(\boldsymbol{U}^{E}_{i,j})\right),0\right\\},$
(32) $\displaystyle
a^{+}_{i+\frac{1}{2},j}=\max\left\\{\lambda_{3K}\left(\frac{\partial\widehat{F}}{\partial\widehat{U}}(\boldsymbol{U}^{W}_{i+1,j})\right),\lambda_{3K}\left(\frac{\partial\widehat{F}}{\partial\widehat{U}}(\boldsymbol{U}^{E}_{i,j})\right),0\right\\},$
$\displaystyle
b^{-}_{i,j+\frac{1}{2}}=\min\left\\{\lambda_{1}\left(\frac{\partial\widehat{G}}{\partial\widehat{U}}(\boldsymbol{U}^{S}_{i,j+1})\right),\lambda_{1}\left(\frac{\partial\widehat{G}}{\partial\widehat{U}}(\boldsymbol{U}^{N}_{i,j})\right),0\right\\},$
$\displaystyle
b^{+}_{i,j+\frac{1}{2}}=\max\left\\{\lambda_{3K}\left(\frac{\partial\widehat{G}}{\partial\widehat{U}}(\boldsymbol{U}^{S}_{i,j+1})\right),\lambda_{3K}\left(\frac{\partial\widehat{G}}{\partial\widehat{U}}(\boldsymbol{U}^{N}_{i,j})\right),0\right\\},$

where $\lambda_{1}$ and $\lambda_{3K}$ are the smallest and the largest
eigenvalues of the matrix, respectively. The semi-discrete version of the
central-upwind scheme is given by

$\displaystyle\dfrac{d}{dt}\overline{\boldsymbol{U}}_{i,j}=-\dfrac{\mathcal{F}_{i+\frac{1}{2},j}-\mathcal{F}_{i-\frac{1}{2},j}}{\Delta
x}-\dfrac{\mathcal{G}_{i,j+\frac{1}{2}}-\mathcal{G}_{i,j-\frac{1}{2}}}{\Delta
y}+\overline{\boldsymbol{S}}_{i,j},$ (33)

where
$\overline{\boldsymbol{S}}_{i,j}\approx\frac{1}{|\mathcal{C}_{i,j}|}\int_{\mathcal{C}_{i,j}}\widehat{S}(\boldsymbol{U},\boldsymbol{B})dxdy$
is the well-balanced discretization of the source term to be discussed in the
following section. The numerical fluxes for $x$\- and $y$-direction are given
by

$\mathcal{F}_{i+\frac{1}{2},j}\coloneqq\dfrac{a^{+}_{i+\frac{1}{2},j}\widehat{F}(\boldsymbol{U}^{E}_{i,j})-a^{-}_{i+\frac{1}{2},j}\widehat{F}(\boldsymbol{U}^{W}_{i+1,j})}{a^{+}_{i+\frac{1}{2},j}-a^{-}_{i+\frac{1}{2},j}}+\dfrac{a^{+}_{i+\frac{1}{2},j}a^{-}_{i+\frac{1}{2},j}}{a^{+}_{i+\frac{1}{2},j}-a^{-}_{i+\frac{1}{2},j}}\left[\boldsymbol{U}^{W}_{i+1,j}-\boldsymbol{U}^{E}_{i,j}\right],$
(34)

and

$\mathcal{G}_{i+\frac{1}{2},j}\coloneqq\dfrac{b^{+}_{i,j+\frac{1}{2}}\widehat{G}(\boldsymbol{U}^{N}_{i,j})-b^{-}_{i,j+\frac{1}{2}}\widehat{G}(\boldsymbol{U}^{S}_{i,j+1})}{b^{+}_{i,j+\frac{1}{2}}-b^{-}_{i,j+\frac{1}{2}}}+\dfrac{b^{+}_{i,j+\frac{1}{2}}b^{-}_{i,j+\frac{1}{2}}}{b^{+}_{i,j+\frac{1}{2}}-b^{-}_{i,j+\frac{1}{2}}}\left[\boldsymbol{U}^{S}_{i,j+1}-\boldsymbol{U}^{N}_{i,j}\right],$
(35)

respectively.

###  4.2 Well-Balanced Flux Discretization

For deterministic 2D SWE, the numerical scheme should accurately capture the
“lake-at-rest” steady state solution, or small perturbations to this state on
coarse meshes. The scheme is said to be well-balanced if it possesses such a
property. A stochastic variant to the “lake-at-rest” property is given by

$q^{x}_{\Lambda}(x,y,t,\xi)=q^{y}_{\Lambda}(x,y,t,\xi)\equiv 0,\quad
h_{\Lambda}+\Pi_{\Lambda}[B](x,y,t,\xi)\equiv C(\xi),$ (36)

where $C(\xi)$ is a random variable depending only on $\xi$. The solution to
(36) represents a still water surface with a flat but stochastic water
surface. The corresponding vector equations for the PCE coefficients are

$\widehat{q^{x}}=\widehat{q^{y}}\equiv\boldsymbol{0},\quad\widehat{h}+\widehat{B}\equiv\widehat{C}.$
(37)

Within the framework of central-upwind scheme, the bottom function
$\widehat{B}$ is first replaced by its continuous piecewise bilinear
interpolant $\widetilde{\widehat{B}}$,

$\displaystyle\widetilde{\widehat{B}}(x,y)$
$\displaystyle=\widehat{B}_{i-\frac{1}{2},j-\frac{1}{2}}+\left(\widehat{B}_{i+\frac{1}{2},j-\frac{1}{2}}-\widehat{B}_{i-\frac{1}{2},j-\frac{1}{2}}\right)\cdot\frac{x-x_{i-\frac{1}{2}}}{\Delta
x}$ $\displaystyle\hskip
55.00008pt+\left(\widehat{B}_{i-\frac{1}{2},j+\frac{1}{2}}-\widehat{B}_{i-\frac{1}{2},j-\frac{1}{2}}\right)\cdot\frac{y-y_{j-\frac{1}{2}}}{\Delta
y}$
$\displaystyle+\left(\widehat{B}_{i+\frac{1}{2},j+\frac{1}{2}}-\widehat{B}_{i+\frac{1}{2},j-\frac{1}{2}}-\widehat{B}_{i-\frac{1}{2},j+\frac{1}{2}}+\widehat{B}_{i-\frac{1}{2},j-\frac{1}{2}}\right)$
$\displaystyle\hskip
55.00008pt\cdot\frac{(x-x_{i-\frac{1}{2}})(y-y_{j-\frac{1}{2}})}{\Delta
x\Delta y},\quad(x,y)\in\mathcal{C}_{i,j}.$

For a general (possibly discontinuous) bottom topography function, the values
at the corners of the cells are estimated by the average of the lower and the
upper limits,

$\displaystyle\widehat{B}_{i\pm\frac{1}{2},j\pm\frac{1}{2}}=\frac{1}{2}\left(\liminf_{(x,y)\to(x_{i\pm\frac{1}{2}},y_{j\pm\frac{1}{2}})}\widehat{B}(x,y)+\limsup_{(x,y)\to(x_{i\pm\frac{1}{2}},y_{j\pm\frac{1}{2}})}\widehat{B}(x,y)\right).$
(38)

For continuous $\widehat{B}(x,y)$, (38) reduces to

$\displaystyle\widehat{B}_{i\pm\frac{1}{2},j\pm\frac{1}{2}}=\widehat{B}(x_{i\pm\frac{1}{2}},y_{j\pm\frac{1}{2}}).$

The cell average of $\widetilde{\widehat{B}}$ in the cell $\mathcal{C}_{i,j}$
is equal to the averages of the values of $\widetilde{\widehat{B}}$ at
midpoints of the edges of the cell, i.e.,

$\displaystyle\overline{\boldsymbol{B}}_{i,j}$
$\displaystyle\coloneqq\frac{1}{|\mathcal{C}_{i,j}|}\iint_{\mathcal{C}_{i,j}}\widetilde{\widehat{B}}(x,y)dxdy,$
$\displaystyle=\frac{1}{4}(\boldsymbol{B}_{i+\frac{1}{2},j}+\boldsymbol{B}_{i-\frac{1}{2},j}+\boldsymbol{B}_{i,j+\frac{1}{2}}+\boldsymbol{B}_{i,j-\frac{1}{2}}),$

where

$\displaystyle\boldsymbol{B}_{i+\frac{1}{2},j}\coloneqq\widetilde{\widehat{B}}(x_{i+\frac{1}{2}},y_{j}),$
$\displaystyle\boldsymbol{B}_{i,j+\frac{1}{2}}\coloneqq\widetilde{\widehat{B}}(x_{i},y_{j+\frac{1}{2}}).$
(39)

Note that, since $\widetilde{\widehat{B}}$ is a piecewise continuous function,
the cell average $\overline{\boldsymbol{B}}_{i,j}$ is equal to the value of
the interpolant at the center of the cell
$\widetilde{\widehat{B}}(x_{i},y_{j})$.

At every time step, we compute the cell averages of the PCE vector of the
water surface
$\overline{\boldsymbol{\eta}}_{i,j}\coloneqq\overline{\boldsymbol{h}}_{i,j}+\overline{\boldsymbol{B}}_{i,j}$
and their second-order-accurate pointwise values
$\boldsymbol{\eta}^{W,E,S,N}_{i,j}$ of the reconstructions. The pointwise
values of the reconstructions for the water heights are then computed by

$\displaystyle\boldsymbol{h}^{W}_{i,j}\coloneqq\boldsymbol{\eta}^{W}_{i,j}-\boldsymbol{B}_{i+\frac{1}{2},j},$
$\displaystyle\boldsymbol{h}^{E}_{i,j}\coloneqq\boldsymbol{\eta}^{E}_{i,j}-\boldsymbol{B}_{i-\frac{1}{2},j},$
(40)
$\displaystyle\boldsymbol{h}^{S}_{i,j}\coloneqq\boldsymbol{\eta}^{S}_{i,j}-\boldsymbol{B}_{i,j-\frac{1}{2}},$
$\displaystyle\boldsymbol{h}^{N}_{i,j}\coloneqq\boldsymbol{\eta}^{N}_{i,j}-\boldsymbol{B}_{i,j+\frac{1}{2}}.$

The numerical fluxes
$\left\\{\left(\mathcal{F}_{i+\frac{1}{2},j},\mathcal{G}_{i,j+\frac{1}{2}}\right)\right\\}_{i\in[N_{x}],j\in[N_{y}]}$
are then computed using the heights defined in (40). With the reconstructions
for the heights and the fluxes, we are now ready to propose the well-balanced
discretization for the source term.

###### Proposition 4.1

The central-upwind scheme (33) is well-balanced if we choose
$\overline{\boldsymbol{S}}_{ij}=\left((\overline{\boldsymbol{S}}^{(1)}_{i,j})^{\top},(\overline{\boldsymbol{S}}^{(2)}_{i,j})^{\top},(\overline{\boldsymbol{S}}^{(3)}_{i,j})^{\top}\right)^{\top}$,
where

$\left\\{\begin{aligned}
\overline{\boldsymbol{S}}^{(1)}_{i,j}&=\boldsymbol{0},\\\
\overline{\boldsymbol{S}}^{(2)}_{i,j}&=-g\mathcal{P}(\overline{\boldsymbol{h}}_{i,j})\left(\dfrac{\boldsymbol{B}_{i+\frac{1}{2},j}-\boldsymbol{B}_{i-\frac{1}{2},j}}{\Delta
x}\right),\\\
\overline{\boldsymbol{S}}^{(3)}_{i,j}&=-g\mathcal{P}(\overline{\boldsymbol{h}}_{i,j})\left(\dfrac{\boldsymbol{B}_{i,j+\frac{1}{2}}-\boldsymbol{B}_{i,j-\frac{1}{2}}}{\Delta
y}\right).\end{aligned}\right.$ (41)

###### Proof 4.1

Let $\boldsymbol{\eta}^{W,E,S,N}_{i,j}$ be the pointwise values of the
reconstructions for the water surface at cell $\mathcal{C}_{i,j}$. Assume at
time $t$, the stochastic water surface is flat and the water is still, that
is,

$\displaystyle\overline{\boldsymbol{\eta}}_{i,j}=\boldsymbol{\eta}^{*}=\text{const
vector},$
$\displaystyle\overline{\boldsymbol{q}^{x}}_{i,j}=\overline{\boldsymbol{q}^{y}}_{i,j}=\boldsymbol{0},$

for all cell $\mathcal{C}_{i,j}$. The piecewise linear reconstructions yield
$\boldsymbol{\eta}^{W,E,S,N}_{i,j}=\boldsymbol{\eta}^{*}$ and
$(\boldsymbol{q}^{x})^{W,E,S,N}_{i,j}=\boldsymbol{0}$. Hence, the numerical
fluxes
$\left\\{\left(\mathcal{F}_{i+\frac{1}{2},j},\mathcal{G}_{i,j+\frac{1}{2}}\right)\right\\}_{i\in[N_{x}],j\in[N_{y}]}$
are

$\mathcal{F}_{i+\frac{1}{2},j}=\begin{pmatrix}\boldsymbol{0}\\\
\frac{g}{2}\mathcal{P}(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i+\frac{1}{2},j})(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i+\frac{1}{2},j})\\\
\boldsymbol{0}\end{pmatrix},$ (42)

and

$\mathcal{G}_{i,j+\frac{1}{2}}=\begin{pmatrix}\boldsymbol{0}\\\
\boldsymbol{0}\\\
\frac{g}{2}\mathcal{P}(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i,j+\frac{1}{2}})(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i,j+\frac{1}{2}})\end{pmatrix}.$
(43)

The semidiscrete form for the stochastic “lake-at-rest” problem is thus given
by

$\displaystyle\dfrac{d}{dt}\overline{\boldsymbol{h}}_{i,j}$
$\displaystyle=\overline{\boldsymbol{S}}^{(1)}_{i,j}\stackrel{{\scriptstyle\eqref{eq:wbproperty}}}{{=}}\boldsymbol{0},$
(44) $\displaystyle\dfrac{d}{dt}\overline{\boldsymbol{q}^{x}}_{i,j}$
$\displaystyle=-\dfrac{1}{\Delta
x}\frac{g}{2}\left[\mathcal{P}(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i+\frac{1}{2},j})(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i+\frac{1}{2},j})-\mathcal{P}(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i-\frac{1}{2},j})(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i-\frac{1}{2},j})\right]+\overline{\boldsymbol{S}}^{(2)}_{i,j},$
$\displaystyle\dfrac{d}{dt}\overline{\boldsymbol{q}^{y}}_{i,j}$
$\displaystyle=-\dfrac{1}{\Delta
y}\frac{g}{2}\left[\mathcal{P}(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i,j+\frac{1}{2}})(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i,j+\frac{1}{2}})-\mathcal{P}(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i,j-\frac{1}{2}})(\boldsymbol{\eta}^{*}-\boldsymbol{B}_{i,j-\frac{1}{2}})\right]+\overline{\boldsymbol{S}}^{(3)}_{i,j},$

The well-balanced property requires that the numerical fluxes to be balanced
by the source term, i.e., the right-hand-side of (44) needs to vanish.
Clearly, $\overline{\boldsymbol{S}}^{(1)}_{i,j}=\boldsymbol{0}$. Using the
linearity of the operator $\mathcal{P}$ and the property (12), one can show
that if we choose $\overline{\boldsymbol{S}}^{(2)}_{i,j}$ and
$\overline{\boldsymbol{S}}^{(3)}_{i,j}$ by (41), the right-hand-side of (44)
will vanish and the well-balanced property of the scheme can be satisfied (see
also [8], Lemma 4.1).

The numerical approximations in (41) can be treated as the approximation of
the cell average of the source term using midpoint quadrature rule with the
partial derivatives
$\boldsymbol{B}_{x}(x_{i},y_{j}),\boldsymbol{B}_{y}(x_{i},y_{j})$ approximated
by the central difference.

###  4.3 Hyperbolicity-Preserving Time Discretization

The hyperbolicity of the system (20) depends only on the vector $\widehat{h}$
and so we focus on the first $K$ equations in (33) in this subsection,

$\displaystyle\dfrac{d}{dt}\overline{\boldsymbol{h}}_{i,j}=-\dfrac{\mathcal{F}^{\widehat{h}}_{i+\frac{1}{2},j}-\mathcal{F}^{\widehat{h}}_{i-\frac{1}{2},j}}{\Delta
x}-\dfrac{\mathcal{G}^{\widehat{h}}_{i,j+\frac{1}{2}}-\mathcal{G}^{\widehat{h}}_{i,j-\frac{1}{2}}}{\Delta
y},$ (45)

where

$\mathcal{F}^{\widehat{h}}_{i+\frac{1}{2},j}=\dfrac{a^{+}_{i+\frac{1}{2},j}(\boldsymbol{q}^{x})_{i,j}^{E}-a^{-}_{i+\frac{1}{2},j}(\boldsymbol{q}^{x})_{i+1,j}^{W}}{a^{+}_{i+\frac{1}{2},j}-a^{-}_{i+\frac{1}{2},j}}+\dfrac{a^{+}_{i+\frac{1}{2},j}a^{-}_{i+\frac{1}{2},j}}{a^{+}_{i+\frac{1}{2},j}-a^{-}_{i+\frac{1}{2},j}}\left[\boldsymbol{h}_{i+1,j}^{W}-\boldsymbol{h}_{i,j}^{E}\right],$
(46)

and

$\mathcal{G}^{\widehat{h}}_{i,j+\frac{1}{2}}=\dfrac{b^{+}_{i,j+\frac{1}{2}}(\boldsymbol{q}^{y})_{i,j}^{N}-b^{-}_{i,j+\frac{1}{2}}(\boldsymbol{q}^{y})_{i,j+1}^{S}}{b^{+}_{i,j+\frac{1}{2}}-b^{-}_{i,j+\frac{1}{2}}}+\dfrac{b^{+}_{i,j+\frac{1}{2}}b^{-}_{i,j+\frac{1}{2}}}{b^{+}_{i,j+\frac{1}{2}}-b^{-}_{i,j+\frac{1}{2}}}\left[\boldsymbol{h}_{i,j+1}^{S}-\boldsymbol{h}_{i,j}^{N}\right].$
(47)

The following lemma presents a computable restriction on a hyperbolicity-
preserving time step for a forward Euler time discretization.

###### Lemma 4.1

Let $\left\\{\xi_{m}\right\\}_{m\in[M]}$ be the nodes of a quadrature rule
satisfying the conditions of 3.2 and
$\boldsymbol{\Phi}(\xi)=(\phi_{1}(\xi),\ldots,\phi_{K}(\xi))^{\top}$. Denote
the numerical approximation to $\overline{\boldsymbol{h}}_{i,j}(t^{n})$ by
$\overline{\boldsymbol{h}}_{i,j}^{n}$ and $\Delta t^{n}\coloneqq
t^{n+1}-t^{n}$. Assume that
$\overline{\boldsymbol{h}}_{i,j}^{n}(\xi_{m})\coloneqq(\overline{\boldsymbol{h}}_{i,j}^{n})^{\top}\Phi(\xi_{m})\;>\;0$
for $m\in[M]$. If $\Delta t^{n}$ satisfies

$\displaystyle\Delta t^{n}$ $\displaystyle<\Delta
t^{n}_{h}\coloneqq\min_{\begin{subarray}{c}m\in[M]\\\
i\in[N_{x}],j\in[N_{y}]\end{subarray}}\left\\{\left|\dfrac{(\overline{\boldsymbol{h}}^{n}_{i,j})^{\top}\boldsymbol{\Phi}(\xi_{m})}{\left[\frac{\mathcal{F}^{\widehat{h}}_{i+\frac{1}{2},j}(t_{n})-\mathcal{F}^{\widehat{h}}_{i-\frac{1}{2},j}(t_{n})}{\Delta
x}+\frac{\mathcal{G}^{\widehat{h}}_{i,j+\frac{1}{2}}(t_{n})-\mathcal{G}^{\widehat{h}}_{i,j-\frac{1}{2}}(t_{n})}{\Delta
y}\right]^{\top}\boldsymbol{\Phi}(\xi_{m})}\right|\right\\},$ (48)

then the flux Jacobian (28),
$J\left(\overline{\boldsymbol{U}}_{i,j}^{n+1}\right)$ is diagonalizable with
real eigenvalues when (20) is discretized using forward Euler solver.

###### Proof 4.2

From 3.2, we only need to guarantee that
$(\overline{\boldsymbol{h}}_{i,j}^{n+1})^{\top}\boldsymbol{\Phi}(\xi_{m})>0$
for all $m\in[M]$, which, for a forward Euler discretization, gives the
inequality

$\displaystyle 0<$
$\displaystyle(\overline{\boldsymbol{h}}^{n+1}_{i,j})^{\top}\boldsymbol{\Phi}(\xi_{m})$
$\displaystyle=$
$\displaystyle(\overline{\boldsymbol{h}}^{n}_{i,j})^{\top}\boldsymbol{\Phi}(\xi_{m})-\frac{\Delta
t^{n}}{\Delta
x}\left[\mathcal{F}^{\widehat{h}}_{i+\frac{1}{2},j}(t_{n})-\mathcal{F}^{\widehat{h}}_{i-\frac{1}{2},j}(t_{n})\right]^{\top}\boldsymbol{\Phi}(\xi_{m})$
(49) $\displaystyle-\frac{\Delta t^{n}}{\Delta
y}\left[\mathcal{G}^{\widehat{h}}_{i,j+\frac{1}{2}}(t_{n})-\mathcal{G}^{\widehat{h}}_{i,j-\frac{1}{2}}(t_{n})\right]^{\top}\boldsymbol{\Phi}(\xi_{m}).$
(50)

Direct calculation and minimization over all the cells as well as all the
quadrature points yields (48).

Combining the hyperbolicity-preserving condition (48) with the wave-speed-
based CFL stability conditions, we choose

$\displaystyle\Delta t^{n}=0.9\min\left\\{\Delta t^{n}_{h},\frac{\min\\{\Delta
x,\Delta y\\}}{2a}\right\\}.$ (51)

Here, $a=\max_{i,j}\\{a_{i+\frac{1}{2},j},b_{i,j+\frac{1}{2}}\\}$ is the
maximum possible magnitude of the propagation speeds at time level $t^{n}$,
where

$\displaystyle
a_{i+\frac{1}{2},j}=\max\\{a^{+}_{i+\frac{1}{2},j},-a^{-}_{i+\frac{1}{2},j}\\},$
$\displaystyle
b_{i,j+\frac{1}{2}}=\max\\{b^{+}_{i,j+\frac{1}{2}},-b^{-}_{i,j+\frac{1}{2}}\\}.$

For now, we assume the local propagation speeds
$a^{\pm}_{i+\frac{1}{2},j},b^{\pm}_{i+\frac{1}{2},j}$ are real, i.e., the SG
system is hyperbolic at the pointwise values of the reconstructions. The
timestep condition analysis for (51) is limited to the cell averages. For
second-order or higher-order central upwind schemes, we need an additional
treatment for the pointwise values $\boldsymbol{U}^{W,E,S,N}_{i,j}$ of the
reconstructions. Further discussion on the procedure to preserve the
hyperbolicity at $\boldsymbol{U}^{W,E,S,N}_{i,j}$ will be provided in Section
4.3.1.

Our numerical scheme uses high-order strong stability-preserving Runge-Kutta
schemes [16] for time discretization, which corresponds to a convex
combination of several forward Euler steps and thus the analysis above for
deriving a hyperbolicity-preserving time step restriction still holds. An
adaptive time control is introduced for the time discretization [5, 19]. Note
that our analysis above can be applied to any other finite volume scheme.

###### Remark 4.1

If the signs of the fluxes at quadrature points are taken into account, the
condition (48) can be relaxed and the scheme will be more efficient.

####  4.3.1 Hyperbolicity-Preserving Filtering for the Pointwise Values of the
Reconstructions

Assuming
$(\overline{\boldsymbol{h}}_{i,j}^{n})^{\top}\boldsymbol{\Phi}(\xi_{m})>0$, we
are able to guarantee that
$(\overline{\boldsymbol{h}}_{i,j}^{n+1})^{\top}\boldsymbol{\Phi}(\xi_{m})>0$
for $m\in[M]$ by the condition (51), provided the local propagation speeds are
real. However, the local propagation speeds are evaluated by the eigenvalues
of the matrices $\frac{\partial\widehat{F}}{\partial\widehat{U}}$ and
$\frac{\partial\widehat{G}}{\partial\widehat{U}}$ (see (23)-(24)) at
$\boldsymbol{U}^{W,E,S,N}_{i,j}$ and, unfortunately, the hyperbolicity at the
pointwise values of the reconstructions is not guaranteed by positivity of
$(\overline{\boldsymbol{h}}_{i,j}^{n})^{\top}\boldsymbol{\Phi}(\xi_{m})$. The
problem can be resolved by filtering the higher-order stochastic moments [37,
8].

Let $p_{\widehat{z}}(\xi)=\sum_{k=1}^{K}\widehat{z}_{k}\phi_{k}(\xi)$ be a
polynomial in $P_{\Lambda}$ with positive moment $\widehat{z}_{1}$. Our goal
is to determine the smallest possible weight $\mu^{\prime}$ such that the
weighted average of the polynomial $p_{\widehat{z}}(\xi)$ and the moment
$\widehat{z}_{1}$ are nonnegative at given quadrature points
$\\{\xi_{m}\\}_{m\in[M]}$, i.e.,

$\mu^{\prime}\widehat{z}_{1}+(1-\mu^{\prime})p_{\widehat{z}}(\xi_{m})\geq
0\Leftrightarrow\widehat{z}_{1}+\sum_{k=2}^{K}(1-\mu^{\prime})\widehat{z}_{k}\phi_{k}(\xi_{m})\geq
0,m\in[M],$ (52)

and we filter the coefficients in the following way,

$\displaystyle\widehat{\mathsf{z}}_{1}=\widehat{z}_{1},$
$\displaystyle\widehat{\mathsf{z}}_{k}=(1-\mu)\widehat{z}_{k},k=2,\ldots,K,$
(53)

where $\mu=\min\\{\mu^{\prime}+\delta,1\\}$, and we select $\delta=10^{-10}$
in our scheme. The filtered polynomial
$p_{\widehat{\mathsf{z}}}(\xi)=\sum_{k=1}^{K}\widehat{\mathsf{z}}_{k}\phi(\xi)$
is positive at quadrature points $\\{\xi_{m}\\}_{m\in[M]}$.

The hyperbolicity at the pointwise values of the reconstructions requires the
positivity of the corresponding reconstructed water height to satisfy
$(\boldsymbol{h}_{i,j}^{W,E,N,S})^{\top}\boldsymbol{\Phi}(\xi_{m})>0$ for all
$m\in[M]$. We will use the filter defined in (52)-(53) to filter the moments
of these polynomials. Our strategy is to take the smallest filtering parameter
$\mu^{n}_{i,j}$, among the filtering parameters of all the pointwise values
$\boldsymbol{h}_{i,j}^{W,E,N,S}$ of the reconstructions in the same cell
$\mathcal{C}_{i,j}$, and use it as the common filtering parameter for
$\boldsymbol{h}_{i,j}^{W,E,N,S}$,

$\displaystyle\left(\mathsf{h}^{W,E,N,S}_{i,j}\right)_{1}=\left(\boldsymbol{h}^{W,E,N,S}_{i,j}\right)_{1},$
(54)
$\displaystyle\left(\mathsf{h}^{W,E,N,S}_{i,j}\right)_{k}=(1-\mu^{n}_{i,j})\left(\boldsymbol{h}^{W,E,N,S}_{i,j}\right)_{k},k=2,\ldots,K.$

The corresponding cell averages are recomputed accordingly,

$\overline{\mathsf{h}}^{n}_{i,j}=\frac{1}{4}\left(\mathsf{h}^{W}_{i,j}+\mathsf{h}^{E}_{i,j}+\mathsf{h}^{N}_{i,j}+\mathsf{h}^{S}_{i,j}\right).$
(55)

There are other choices of filters to further reduce the oscillations (e.g.,
see [26, 38, 27]), but we adopt this filter because of its simplicity. Note
that, in (54)-(55), we only filter the polynomial moments for
$\boldsymbol{h}^{W,E,N,S}$ and do not filter the discharges.

The above strategy requires that the first moments of
$\boldsymbol{h}_{i,j}^{W,E,N,S}$ are positive during computation, but this is
not necessarily the case. This may lead to the failure of the filtering
procedure (54). This is an issue that happens even for the deterministic
system of 2D SWE [23]. Therefore, we propose the following correction for the
PCE vector of the reconstructed water height with nonpositive first moment,

$\displaystyle\textrm{if }\left(\boldsymbol{h}^{W}_{i,j}\right)_{1}\leq
0\hskip 5.0pt\textrm{then take
}\boldsymbol{h}^{W}_{i,j}=\boldsymbol{0},\;\;\boldsymbol{h}^{E}_{i,j}=2\overline{\boldsymbol{h}}^{n}_{i,j},$
$\displaystyle\textrm{if }\left(\boldsymbol{h}^{E}_{i,j}\right)_{1}\leq
0\hskip 5.0pt\textrm{then take
}\boldsymbol{h}^{E}_{i,j}=\boldsymbol{0},\;\;\boldsymbol{h}^{W}_{i,j}=2\overline{\boldsymbol{h}}^{n}_{i,j},$
$\displaystyle\textrm{if }\left(\boldsymbol{h}^{N}_{i,j}\right)_{1}\leq
0\hskip 5.0pt\textrm{then take
}\boldsymbol{h}^{N}_{i,j}=\boldsymbol{0},\;\;\boldsymbol{h}^{S}_{i,j}=2\overline{\boldsymbol{h}}^{n}_{i,j},$
$\displaystyle\textrm{if }\left(\boldsymbol{h}^{S}_{i,j}\right)_{1}\leq
0\hskip 5.0pt\textrm{then take
}\boldsymbol{h}^{N}_{i,j}=\boldsymbol{0},\;\;\boldsymbol{h}^{N}_{i,j}=2\overline{\boldsymbol{h}}^{n}_{i,j}.$

In this way, the cell averages remain consistent, i.e., the cell average is
equal to the average of the pointwise values for the corresponding cell. Note
that, in the case when the first moment is $0$, the filter will not be
applied. This correction reduces to a similar procedure in the central-upwind
scheme for the deterministic shallow water equations [23].

####  4.3.2 Velocity Desingularization

In computation, the matrix $\mathcal{P}(\boldsymbol{h}^{W,E,N,S}_{i,j})$ may
be near singular at the pointwise values of the reconstructions and may result
in large round-off errors when determining the PCE vector for the velocity
(25). An approach for stochastic velocity desingularization is proposed in
[8]. To recall the idea, we consider $\mathcal{P}(\boldsymbol{h}^{W}_{i,j})$
as an example. But our process can be extended naturally to the other point
values of the reconstructions. Let

$\displaystyle\mathcal{P}(\boldsymbol{h}^{W}_{i,j})=Q^{\top}\Sigma Q,$ (56)

be the eigenvalues decomposition of $\mathcal{P}(\boldsymbol{h}^{W}_{i,j})$,
where $\Sigma=\text{diag}(\sigma_{1},\ldots,\sigma_{K})$. For $k\in[K]$ and a
prescribed $\epsilon$, defined the corrected inverse of $\Sigma$ to be,

$\displaystyle\Sigma^{\text{cor}}$
$\displaystyle=\text{diag}(\sigma^{\text{cor}}_{1},\ldots,\sigma^{\text{cor}}_{K}),$
$\displaystyle\sigma^{\text{cor}}_{k}=\frac{\sqrt{2}\sigma_{k}}{\sqrt{\sigma^{4}_{k}+\max\\{\sigma^{4}_{k},\epsilon^{4}\\}}}.$
(57)

In our simulation, $\epsilon=\min\\{\Delta x,\Delta y\\}$. The PCE
coefficients for the velocity vector are then computed by,

$\displaystyle\boldsymbol{u}^{W}_{i,j}=Q^{\top}\Sigma^{\text{cor}}Q(\boldsymbol{q}^{x})^{W}_{i,j},$
$\displaystyle\boldsymbol{v}^{W}_{i,j}=Q^{\top}\Sigma^{\text{cor}}Q(\boldsymbol{q}^{y})^{W}_{i,j}.$
(58)

For well-conditioned $\mathcal{P}(\boldsymbol{h}^{W}_{i,j})$, the correction
(58) reduces to (25) and no further reconstruction is needed. In the case when
$\mathcal{P}(\boldsymbol{h}^{W}_{i,j})$ is ill-conditioned, to remain the
consistency of the PCE vector of the discharges, we recompute the discharges
using the corrected velocity

$\displaystyle(\boldsymbol{q}^{x})^{W}_{i,j}=\mathcal{P}(\boldsymbol{h}^{W}_{i,j})\boldsymbol{u}^{W}_{i,j},$
$\displaystyle(\boldsymbol{q}^{y})^{W}_{i,j}=\mathcal{P}(\boldsymbol{h}^{W}_{i,j})\boldsymbol{v}^{W}_{i,j}.$
(59)

If there is no uncertainty in the system (all moments beyond the first
vanish), the desingularization procedure reduces to the deterministic velocity
desingularization procedure proposed in [23, 19]. Readers can also find other
desingularization options there.

##  5 Numerical Results

In this section, we will demonstrate the capability of our scheme to capture
the solutions to small perturbations of the stochastic “lake-at-rest” solution
for different choices of the random variables $\xi$. The distributions we
consider are from the families of distributions with Beta density. For a
single random variable $\xi^{(1)}$, its density is given by,

$\displaystyle\rho(\xi^{(1)})\coloneqq\rho^{(\alpha,\beta)}(\xi^{(1)})$
$\displaystyle=C(\alpha,\beta)(1-\xi^{(1)})^{\alpha}(1+\xi^{(1)})^{\beta},$
$\displaystyle C(\alpha,\beta)^{-1}$
$\displaystyle=2^{\alpha+\beta+1}\mathcal{B}(\beta+1,\alpha+1),$ (60)

where $\mathcal{B}(\cdot,\cdot)$ is the Beta function, and the parameters
$\alpha,\beta>-1$ can be chosen freely and control how mass concentrates at
$\xi^{(1)}=1$ and $\xi^{(1)}=-1$, respectively. The associated orthonormal
polynomial basis is the Jacobi polynomials with parameter $\alpha$ and
$\beta$. In particular, $\alpha=\beta=0$ corresponds to the uniform
distribution $\mathcal{U}(-1,1)$ and the Legendre polynomial basis.

In our numerical tests, we consider the cases with both the one-dimensional
random variable $\xi=\xi^{(1)}$ and the two-dimensional random variable
$\xi=(\xi^{(1)},\xi^{(2)})^{\top}$, where $\xi^{(1)}$ and $\xi^{(2)}$ are
independent one-dimensional random variables with different distributions, and
thus the density function for $\xi$ is defined by
$\rho(\xi)\coloneqq\rho(\xi^{(1)})\rho(\xi^{(2)})$. The multi-index set
$\Lambda$ associated with the polynomial space $P_{\Lambda}$ is chosen to be
the tensor-product set,

$\displaystyle\Lambda=\\{(\nu^{(1)},\nu^{(2)})\in\mathbbm{N}^{2}\;\;|\;\;0\leq\nu^{(1)},\nu^{(2)}\leq
3\\},$ (61)

for all the two-dimensional random field numerical examples.

Throughout the numerical section, the gravitational constant is $g=1$, the
parameter for the generalized minmod limiter is set to be $\theta=1.3$ (except
for Section 5.4, where $\theta=1.0$). The filtering parameter is set to be
$\delta=10^{-10}$.

###  5.1 Example 1: Accuracy Test

This subsection is devoted to checking the order of accuracy of the proposed
scheme for our SG system of 2D SWE. We consider a stochastic version of the
accuracy test in [4] with deterministic initial water surface,

$w(x,y,0,\xi)=1,\quad u(x,y,0,\xi)=0.3,\quad v(x,y,0,\xi)=0,$ (62)

and stochastic elliptic-shaped hump bottom,

$B(x,y)=0.5e^{-25(x-1)^{2}-50(y-0.5)^{2}}+0.1(\xi+1).$ (63)

The one-dimensional random variable is chosen to be
$\xi\sim\mathcal{U}(-1,1)$. The computational domain is $[0,2]\times[0,1]$.
All the boundary conditions are obtained by zero-order extrapolation. The
reference solution is computed on a $800\times 800$ uniform rectangular grid.
By $t=0.07$, the solution converges to the steady state, which, for this
example, is non constant but smooth. The contour plots of the reference
solution for the mean water surface and the standard deviation of the water
surface at $t=0.07$ are shown in Figure 2.

(a)

(b)

Figure 2: Contours for the water surface of the reference solution. Left:
mean. Right: standard deviation.

We computed the error between the reference solution and the numerical
solution using,

$\displaystyle\text{err}(U_{d})=\|U_{d}(x,y,\xi)-U_{\text{ref}}(x,y,\xi)\|_{L^{1}(\Omega_{xy};L^{2}_{\rho}(\Omega_{\xi};\mathbbm{R}^{3}))},$
(64)

where $\Omega_{xy}$ is the physical domain, $\Omega_{\xi}$ is the stochastic
domain, $U_{d}$ is the numerical $K$-term PCE solution to the
$\xi$-parameterized equation (1), and $U_{\text{ref}}$ is the $K$-term PCE
reference solution to the $\xi$-parameterized equation (1). Here, for the
stochastic state variable in $K$-term PCE,
$U=(h_{\Lambda},q^{x}_{\Lambda},q^{y}_{\Lambda})^{\top}$, the norm
$\|\cdot\|_{L^{1}(\Omega_{xy};L^{2}_{\rho}(\Omega_{\xi};\mathbbm{R}^{3}))}$ is
defined to be,

$\displaystyle\|U(x,y,t,\xi)\|_{L^{1}(\Omega_{xy};L^{2}_{\rho}(\Omega_{\xi};\mathbbm{R}^{3}))}\coloneqq$
$\displaystyle\int_{\Omega_{xy}}\left|\|h_{\Lambda}(x,y,t,\xi)\|_{L^{2}_{\rho}}+\|q^{x}_{\Lambda}(x,y,t,\xi)\|_{L^{2}_{\rho}}+\|q^{y}_{\Lambda}(x,y,t,\xi)\|_{L^{2}_{\rho}}\right|dx\;dy.$

From the orthogonality of the basis functions in the PCE, we have,

$\|h_{\Lambda}(x,y;\xi)\|_{L^{2}_{\rho}}=\|\widehat{h}\|_{2},$

where $\widehat{h}$ is the vector of the PCE coefficients of $h_{\Lambda}$,
and $\|\cdot\|_{2}$ is the Euclidean norm. The computations for the norm of
$q^{x}_{\Lambda}$ and $q^{y}_{\Lambda}$ are similar. The orders of convergence
for the numerical solutions with $K=4$ and $K=8$ are reported in Table 1. The
orders are similar to the orders of the deterministic accuracy test reported
in [4] and we expect that using a finer grid for the reference solution will
yield second-order convergence. We also note that, for this example,
increasing the polynomial order will not have a significant effect on the
order of the convergence, especially when $K$ is relatively large ($K>8$).
Most high order moments remain to be near $0$ at $t=0.07$.

Grid size | K=4 | K=8
---|---|---
Error | Order | Error | Order
$100\times 100$ | 1.475875e-05 | $-$ | 1.475865e-05 | $-$
$200\times 200$ | 4.343711e-06 | 1.764571 | 4.343714e-06 | 1.764559
$400\times 400$ | 1.296122e-06 | 1.744727 | 1.296122e-06 | 1.744727
Table 1: Order of convergence for (62)-(63) with different grid sizes and
different $K$ values.

###  5.2 Example 2: Comparison with Collocation Method

The second numerical example serves to compare our SG formulations (20) as
well as our scheme with a collocation procedure. Consider a deterministic
initial water surface

$w(x,y,0,\xi)=\left\\{\begin{aligned} &1.01,&&\text{if }0.05<x<0.15,\\\
&1,&&\text{otherwise},\end{aligned}\right.\quad u(x,y,0,\xi)=v(x,y,0,\xi)=0,$
(65)

with an uncertain bottom function

$B(x,y,\xi)=0.8e^{-5(x-0.9)^{2}-50(y-0.5)^{2}}+0.1(\xi+1),$ (66)

where $\xi$ has Beta density with parameters $(\alpha,\beta)=(1,3)$. At $t=0$,
the water surface is flat except at $0.05<x<0.15$, and the uncertainty of the
bottom skews toward the positive direction, i.e., the probability for $\xi$ to
be a positive number is higher. The computational domain is
$[0,2]\times[0,1]$. The boundary conditions for the upper and the lower
boundaries are periodic boundary conditions. The boundary conditions for the
left and right boundaries are outflow boundary conditions. We compute the SG
solution for $K=4,8$ on $200\times 200$ and $400\times 400$ grid up to time
$T=1.2$. As a comparison, we compute the stochastic collocation (SC) solution
$U_{SC}$ for $K=8$ with $M=100$ collocation points on $400\times 400$ grid for
the same terminal time, where the SC solution is obtained by

$\displaystyle
U_{SC}(x,y,t,\zeta)=\sum_{k=1}^{K}\widehat{U}_{SC,k}(x,y,t)\phi_{k}(\zeta),$
$\displaystyle\widehat{U}_{SC,k}(x,y,t)\coloneqq\sum_{m=1}^{M}U_{SC}(x,y,t,\zeta_{m})\phi_{k}(\zeta_{m})z_{m}.$

Here, $\\{\zeta_{m},z_{m}\\}_{m\in[M]}$ is an $M$-point $\rho$-Gaussian
quadrature rule. The comparison of the contours of the statistics of the water
surface is shown in Figure 3 (mean) and Figure 4 (standard deviation). Since
the SG solutions and the SC solutions are obtained from numerically solving
different partial differential equations, we do not expect that one solution
will converge to the other solution by, for example, increasing the number $M$
of quadrature points and/or refining the mesh. The contour lines look
different near the intermediate region from peak to valley of the mean water
surface, but overall they are similar, especially between the SG solutions on
$400\times 400$ grid with $K=8$ and the SC solution.

(a)

(b)

(c)

(d)

(e)

Figure 3: Numerical solutions to (65)-(66), mean water surface at $t=1.2$,
grid size $400\times 400$. Top left: SG solution, $K=4$, $200\times 200$ grid;
Top right: SG solution, $K=8$, $200\times 200$ grid. Middle left: SG solution,
$K=4$, $400\times 400$ grid; Middle right: SG solution, $K=8$, $400\times 400$
grid. Bottom: SC solution, $K=8$, $400\times 400$ grid. The number of
collocation points is $100$.

(a)

(b)

(c)

(d)

(e)

Figure 4: Numerical solutions to (65)-(66), standard deviation of the water
surface at $t=1.2$.Top left: SG solution, $K=4$, $200\times 200$ grid; Top
right: SG solution, $K=8$, $200\times 200$ grid. Middle left: SG solution,
$K=4$, $400\times 400$ grid; Middle right: SG solution, $K=8$, $400\times 400$
grid. Bottom: SC solution, $K=8$, $400\times 400$ grid. The number of
collocation points is $100$.

###  5.3 Example 3: A Two-Dimensional Random Variable Example

In this example, we will consider a more difficult variant of (65)-(66).
Consider a deterministic initial water surface

$w(x,y,0,\xi)=\left\\{\begin{aligned} &1.01,&&\text{if }0.05<x<0.15,\\\
&1,&&\text{otherwise},\end{aligned}\right.\quad u(x,y,0,\xi)=v(x,y,0,\xi)=0,$
(67)

and a stochastic bottom

$B(x,y,\xi)=0.8e^{-5(x-0.9+0.1\xi^{(1)})^{2}-50(y-0.5+0.1\xi^{(2)})^{2}},$
(68)

where $\xi=(\xi^{(1)},\xi^{(2)})^{\top}$ is a two-dimensional random variable,
with $\xi^{(1)}$ having Beta density with parameters $(\alpha,\beta)=(1,3)$,
and $\xi^{(2)}\sim\mathcal{U}(-1,1)$. The randomness is put on the position of
the center of the Gaussian-shape “hump”.

The computational domain is $[0,2]\times[0,1]$. We compute the numerical
solution on $100\times 100$ and $200\times 200$ grid at various time. We plot
the contours of mean overlaid with disk glyphs (Figure 5 and Figure 6) whose
radius are proportional to the uncertainties at the corresponding cells. We
uniformly downsample the cells for the disk glyphs to better present the
uncertainties in the figures. The water propagates towards the right in the
computational domain and the uncertainties remain small for most of the
sampled cells at the beginning ($t=0.6$). As the water interacts with the
hump, the water splits and propagates along all directions. The contours of
the mean water surface generate more complicated structures and the
uncertainties propagate along the directions of the mean waves. The largest
uncertainties were concentrated near the peak of the water surface
($t=0.6,0.9$) but then gradually “spread out” following the mean waves
($t=1.2,1.5,1.8$). One can also compare the middle left contours
(corresponding to the solution at $t=1.2$) in Figure 6 or Figure 5 with the
contours of the mean water surfaces in Figure 3. The uncertainties imposed in
the position shifts the peak of the mean towards the right. It is reasonable
because $\xi^{(1)}$ follows a distribution that skews negatively so the mean
center of the hump $(0.9-0.1\xi^{(1)})$ is skews towards the right. On the
other hand, the uncertainties on the center increase the width of the peak and
valley structures.

(a)

(b)

(c)

(d)

(e)

Figure 5: Numerical solution to (67)-(68), water surface, $100\times 100$
grid, disk-glyph over mean contours, where the radii of the disks indicate the
magnitude of the standard deviation, $t=0.6,0.9,1.2,1.5,1.8$ (Top left, top
right, middle left, middle right, bottom). The largest disks are corresponding
to the standard deviation values 6.97e-4, 1.50e-3, 1.03e-3, 5.16e-4, and
4.81e-4, respectively.

(a)

(b)

(c)

(d)

(e)

Figure 6: Numerical solution to (67)-(68), $200\times 200$ grid, water
surface, disk-glyph over mean contours, where the radii of the disks indicate
the magnitude of the standard deviation, $t=0.6,0.9,1.2,1.5,1.8$ (Top left,
top right, middle left, middle right, bottom). The largest disks are
corresponding to the standard deviation values 2.20e-3, 3.35e-3, 2.00e-3,
1.58e-3, and 1.20e-3, respectively.

###  5.4 Example 4: A Submerged Flat Plateau

We consider a deterministic initial water surface with deterministic
perturbation

$w(x,y,0)=\left\\{\begin{aligned} &1.0001,&&\text{if }-0.4<x<-0.3,\\\
&1,&&\text{otherwise},\end{aligned}\right.,\quad u(x,y,0)=v(x,y,0)=0,$ (69)

and with stochastic bottom topography,

$B(x,y,\xi)=\left\\{\begin{aligned} &0.9998,&&r\leq 0.1,\\\
&9.997(0.2-r)+0.0001(\xi^{(1)}+1),&&0.1<r\leq 0.2,\\\
&0.0001,&&\text{otherwise},\end{aligned}\right.$ (70)

where $r\coloneqq\sqrt{x^{2}+y^{2}}+0.0001(\xi^{(2)}+1)$. The bottom is a flat
plateau beneath, but very close to, the water surface. It is continuous but
not smooth in the physical domain. The uncertainties
$\xi^{(1)}\sim\mathcal{U}(-1,1)$ are on the heights of the bottom in the
transition region between the plateau and the lowland in the water. The
uncertainties $\xi^{(2)}$ (with Beta density $(\alpha,\beta)=(1,3)$) are on
the radius of the disk-shaped plateau. We compute the solutions at various
time on the computational domain $[-0.5,0.5]\times[-0.5,0.5]$ with a
$100\times 100$ and $200\times 200$ grid. The mean contour plots with disk
glyphs are presented in Figure 7 and Figure 8. The initial deterministic
perturbation is split into two waves. The left-going wave propagates out of
the computational region while the right-going wave interacts with the plateau
and generates complicated waves as the interaction starts. The initial
deterministic perturbation on the water surface is $10^{-4}$ and largest level
of the mean water displacement from the steady state stay in the same order
during the computation. The uncertainties (standard deviation) remain small
and are close to the island, which shows the capability of our scheme to
capture the stochastic “lake-at-rest” solution (36).

(a)

(b)

(c)

(d)

Figure 7: Numerical solution to (69)-(70), $100\times 100$, water surface,
disk-glyph over mean contours, where the radii of the disks indicate the
magnitude of the standard deviation, $t=0.2,0.35,0.5,0.65$ (top left, top
right, middle left, middle right). The largest disks are corresponding to the
standard deviation values 1.27e-8, 5.32e-7, 3.60e-7, and 5.37e-7,
respectively.

(a)

(b)

(c)

(d)

Figure 8: Numerical solution to (69)-(70), $200\times 200$, water surface,
disk-glyph over mean contours, where the radii of the disks indicate the
magnitude of the standard deviation, $t=0.2,0.35,0.5,0.65$ (top left, top
right, middle left, middle right). The largest disks are corresponding to the
standard deviation values 4.83e-8, 6.57e-6, 4.20e-6, and 3.50e-6,
respectively.

###  5.5 Example 5: Another Two-Dimensional Example

We consider the deterministic initial water surface

$w(x,y,0,\xi)=1,\quad u(x,y,0,\xi)=0.3,\quad v(x,y,0,\xi)=0,$ (71)

with stochastic bottom

$B(x,y,\xi)=0.5e^{-12.5(\xi^{(1)}+1)(x-1)^{2}-25(\xi^{(2)}+1)(y-0.5)^{2}}.$
(72)

where $\xi=(\xi^{(1)},\xi^{(2)})$ is a random vector, $\xi^{(1)}$ is the
distribution with Beta density $(\alpha,\beta)=(3,1)$,
$\xi^{(2)}\sim\mathcal{U}(-1,1)$. This example is a multidimensional variant
of Section 5.1. Here, we put the uncertainties on the width of the Gaussian-
shape hump. The computational domain is $[0,2]\times[0,1]$. We compute the
solution numerically on $100\times 100$ (Figure 9) and $200\times 200$ (Figure
10) grids up to $t=0.07$. Compared with the solution in Figure 2, where only a
simple additive uncertainty is imposed, the contours in Figure 9 (or Figure
10) are more complicated with the introduction of the two dimensional
uncertainties on the width of the hump. The largest uncertainties occur near
the peak and the lowest part of the mean waves. The region between them,
however, has very small uncertainties.

(a)

(b)

(c)

Figure 9: Numerical solution to (71)-(72), water surface $t=0.07$, $100\times
100$ grid. Top: mean; Bottom: disk-glyph over mean contours, where the radii
of the disks indicate the magnitude of the standard deviation. The largest
disk is corresponding to the standard deviation 8.10e-3.

(a)

(b)

(c)

Figure 10: Numerical solution to (71)-(72), water surface $t=0.07$, $200\times
200$ grid. Top: mean; Bottom: disk-glyph over mean contours, where the radii
of the disks indicate the magnitude of the standard deviation. The largest
disk is corresponding to the standard deviation 8.58e-3.

##  6 Conclusion

We have introduced a hyperbolicity-preserving SG formulations for the two-
dimensional shallow water equations. A sufficient condition for preserving the
hyperbolicity of our SG system has been introduced and can be interpreted as a
“stochastic” positivity condition. Our formulation remains in conservative
form. Based on our theoretical results, we develop a well-balanced
hyperbolicity-preserving central-upwind solver and illustrate the
effectiveness of our solver on some challenging numerical examples. However,
our theoretical framework is not limited to central-upwind type schemes, and
can be extended to other finite-volume-type schemes (e.g., [39]) or
discontinuous Galerkin schemes ([46]) with modifications on the discretization
of the well-balanced source term, as well as the CFL-type time condition.

Currently, our scheme is not designed to handle the dry/wet interfaces, for
example, a wave arriving or leaving the shore. The presence of dry or
partially dry cells may cause computationally impractical time step even for
deterministic problems. In deterministic cases, the “draining” time step
strategy has been introduced to capture the wet/dry interfaces [3, 2, 30].
Future design of stochastic version of “draining” time step may adapt our
scheme to more challenging problems. Another possible extension is to preserve
the strict positivity of the water height rather than the positivity at a
finite number of quadrature points, which requires an efficient algorithm to
preserve the positivity of high-order polynomials (e.g., [51]).

## Acknowledgement

A. Narayan was partially supported by NSF DMS-1848508.

## References

  * Babuska et al. [2004] Babuska, I., Tempone, R., Zouraris, G.E., 2004.  Galerkin Finite Element Approximations of Stochastic Elliptic Partial Differential Equations.  SIAM Journal on Numerical Analysis 42, 800–825.  doi:https://doi.org/10.1137/S0036142902418680.
  * Bollermann et al. [2013] Bollermann, A., Chen, G., Kurganov, A., Noelle, S., 2013\.  A Well-Balanced Reconstruction of Wet/Dry Fronts for the Shallow Water Equations.  Journal of Scientific Computing 56, 267–290.  doi:https://doi.org/10.1007/s10915-012-9677-5.
  * Bollermann et al. [2011] Bollermann, A., Noelle, S., Lukáčová-Medvid’ová, M., 2011\.  Finite Volume Evolution Galerkin Methods for the Shallow Water Equations with Dry Beds.  Communications in Computational Physics 10, 371–404.  doi:https://doi.org/10.4208/cicp.220210.020710a.
  * Bryson et al. [2011] Bryson, S., Epshteyn, Y., Kurganov, A., Petrova, G., 2011\.  Well-balanced positivity preserving central-upwind scheme on triangular grids for the saint-venant system.  ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique 45, 423–446.  doi:https://doi.org/10.1051/m2an/2010060.
  * Chertock et al. [2015a] Chertock, A., Cui, S., Kurganov, A., Wu, T., 2015a.  Well-balanced positivity preserving central-upwind scheme for the shallow water system with friction terms.  International Journal for Numerical Methods in Fluids 78, 355–383.  doi:https://doi.org/10.1002/fld.4023.
  * Chertock et al. [2015b] Chertock, A., Jin, S., Kurganov, A., 2015b.  A well-balanced operator splitting based stochastic Galerkin method for the one-dimensional Saint-Venant system with uncertainty.  preprint URL: https://ki-net.umd.edu/pubs/files/SW.pdf.
  * Chertock et al. [2015c] Chertock, A., Jin, S., Kurganov, A., 2015c.  An operator splitting based stochastic Galerkin method for the one-dimensional compressible Euler equations with uncertainty.  preprint , 1–21URL: http://www.math.wisc.edu/~jin/PS/Euler-UQ.pdf.
  * Dai et al. [2021] Dai, D., Epshteyn, Y., Narayan, A., 2021.  Hyperbolicity-Preserving and Well-Balanced Stochastic Galerkin Method for Shallow Water Equations.  SIAM Journal on Scientific Computing 43, A929–A952.  doi:https://doi.org/10.1137/20M1360736.
  * Debusschere et al. [2004] Debusschere, B.J., Najm, H.N., Pébay, P.P., Knio, O.M., Ghanem, R.G., Le Maître, O.P., 2004\.  Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes.  SIAM journal on scientific computing 26, 698–719.  doi:https://doi.org/10.1137/S1064827503427741.
  * Després et al. [2013] Després, B., Poëtte, G., Lucor, D., 2013.  Robust Uncertainty Propagation in Systems of Conservation Laws with the Entropy Closure Method, in: Uncertainty Quantification in Computational Fluid Dynamics. Springer, pp. 105–149.
  * Eigel et al. [2014] Eigel, M., Gittelson, C.J., Schwab, C., Zander, E., 2014\.  Adaptive stochastic Galerkin FEM.  Computer Methods in Applied Mechanics and Engineering 270, 247–269.  doi:https://doi.org/10.1016/j.cma.2013.11.015.
  * Ernst et al. [2012] Ernst, O.G., Mugler, A., Starkloff, H.J., Ullmann, E., 2012\.  On the Convergence of Generalized Polynomial Chaos Expansions.  ESAIM: Mathematical Modelling and Numerical Analysis 46, 317–339.  doi:https://doi.org/10.1051/m2an/2011045.
  * Gerster and Herty [2020] Gerster, S., Herty, M., 2020\.  Entropies and Symmetrization of Hyperbolic Stochastic Galerkin Formulations.  Communications in Computational Physics 27, 639–671.  doi:https://doi.org/10.4208/cicp.OA-2019-0047.
  * Gerster et al. [2019] Gerster, S., Herty, M., Sikstel, A., 2019.  Hyperbolic stochastic Galerkin formulation for the p-system.  Journal of Computational Physics 395, 186–204.  doi:https://doi.org/10.1016/j.jcp.2019.05.049.
  * Ghanem and Spanos [1991] Ghanem, R.G., Spanos, P.D., 1991\.  Stochastic Finite Elements: A Spectral Approach .
  * Gottlieb et al. [2001] Gottlieb, S., Shu, C.W., Tadmor, E., 2001.  Strong Stability-Preserving High-Order Time Discretization Methods.  SIAM review 43, 89–112.  doi:https://doi.org/10.1137/S003614450036757X.
  * Hu and Jin [2016] Hu, J., Jin, S., 2016.  A stochastic Galerkin method for the Boltzmann equation with uncertainty.  Journal of Computational Physics 315, 150–168.  doi:https://doi.org/10.1016/j.jcp.2016.03.047.
  * Jin and Shu [2019] Jin, S., Shu, R., 2019.  A Study of Hyperbolicity of Kinetic Stochastic Galerkin System for the Isentropic Euler Equations with Uncertainty.  Chinese Annals of Mathematics, Series B 40, 765–780.  doi:https://doi.org/10.1007/s11401-019-0159-z.
  * Kurganov [2018] Kurganov, A., 2018.  Finite-volume schemes for shallow-water equations.  Acta Numerica 27, 289–351.  doi:https://doi.org/10.1017/S0962492918000028.
  * Kurganov and Levy [2002] Kurganov, A., Levy, D., 2002\.  Central-upwind schemes for the Saint-Venant system.  ESAIM: Mathematical Modelling and Numerical Analysis 36, 397–425.  doi:https://doi.org/10.1051/m2an:2002019.
  * Kurganov and Lin [2007] Kurganov, A., Lin, C.T., 2007\.  On the Reduction of Numerical Dissipation in Central-Upwind Schemes.  Communications in Computational Physics 2, 141–163.
  * Kurganov et al. [2001] Kurganov, A., Noelle, S., Petrova, G., 2001.  Semidiscrete Central-Upwind Schemes for Hyperbolic Conservation Laws and Hamilton–Jacobi Equations.  SIAM Journal on Scientific Computing 23, 707–740.  doi:https://doi.org/10.1137/S1064827500373413.
  * Kurganov and Petrova [2007] Kurganov, A., Petrova, G., 2007\.  A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system.  Communications in Mathematical Sciences 5, 133–160.
  * Kurganov et al. [2007] Kurganov, A., Petrova, G., Popov, B., 2007.  Adaptive Semidiscrete Central-Upwind Schemes for Nonconvex Hyperbolic Conservation Laws.  SIAM Journal on Scientific Computing 29, 2381–2401.  doi:https://doi.org/10.1137/040614189.
  * Kurganov and Tadmor [2000] Kurganov, A., Tadmor, E., 2000\.  New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations.  Journal of Computational Physics 160, 241–282.  doi:https://doi.org/10.1006/jcph.2000.6459.
  * Kusch et al. [2020] Kusch, J., McClarren, R.G., Frank, M., 2020.  Filtered stochastic Galerkin methods for hyperbolic equations.  Journal of Computational Physics 403, 109073.  doi:https://doi.org/10.1016/j.jcp.2019.109073.
  * Kusch and Schlachter [2020] Kusch, J., Schlachter, L., 2020\.  Oscillation Mitigation of Hyperbolicity-Preserving Intrusive Uncertainty Quantification Methods for Systems of Conservation Laws.  arXiv:2008.07845 .
  * Le Maître et al. [2004] Le Maître, O., Knio, O., Najm, H., Ghanem, R., 2004\.  Uncertainty propagation using Wiener–Haar expansions.  Journal of computational Physics 197, 28–57.  doi:https://doi.org/10.1016/j.jcp.2003.11.033.
  * Le Maître and Knio [2010] Le Maître, O., Knio, O.M., 2010\.  Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics.  Springer Science & Business Media.
  * Liu et al. [2018] Liu, X., Albright, J., Epshteyn, Y., Kurganov, A., 2018\.  Well-balanced positivity preserving central-upwind scheme with a novel wet/dry reconstruction on triangular grids for the Saint-Venant system.  Journal of Computational Physics 374, 213–236.  doi:https://doi.org/10.1016/j.jcp.2018.07.038.
  * Mishra et al. [2012] Mishra, S., Schwab, C., Sukys, J., 2012.  Multilevel Monte Carlo Finite Volume Methods for Shallow Water Equations with Uncertain Topography in Multi-dimensions.  SIAM Journal on Scientific Computing 34, B761–B784.  doi:https://doi.org/10.1137/110857295.
  * Nessyahu and Tadmor [1990] Nessyahu, H., Tadmor, E., 1990\.  Non-oscillatory central differencing for hyperbolic conservation laws.  Journal of computational physics 87, 408–463.  doi:https://doi.org/10.1016/0021-9991(90)90260-8.
  * Nobile et al. [2008] Nobile, F., Tempone, R., Webster, C.G., 2008.  A Sparse Grid Stochastic Collocation Method for Partial Differential Equations with Random Input Data.  SIAM Journal on Numerical Analysis 46, 2309–2345.  doi:https://doi.org/10.1137/060663660.
  * Pettersson et al. [2014] Pettersson, P., Iaccarino, G., Nordström, J., 2014.  A stochastic Galerkin method for the Euler equations with Roe variable transformation.  Journal of Computational Physics 257, 481–500.  doi:https://doi.org/10.1016/j.jcp.2013.10.011.
  * Poëtte [2019] Poëtte, G., 2019.  Contribution to the mathematical and numerical analysis of uncertain systems of conservation laws and of the linear and nonlinear Boltzmann equation.  Ph.D. thesis.  URL: https://hal.archives-ouvertes.fr/tel-02288678/document.
  * Poëtte et al. [2009] Poëtte, G., Després, B., Lucor, D., 2009.  Uncertainty quantification for systems of conservation laws.  Journal of Computational Physics 228, 2443–2467.  doi:https://doi.org/10.1016/j.jcp.2008.12.018.
  * Schlachter and Schneider [2018] Schlachter, L., Schneider, F., 2018\.  A hyperbolicity-preserving stochastic Galerkin approximation for uncertain hyperbolic systems of equations.  Journal of Computational Physics 375, 80–98.  doi:https://doi.org/10.1016/j.jcp.2018.07.026.
  * Schlachter et al. [2020] Schlachter, L., Schneider, F., Kolb, O., 2020.  Weighted Essentially Non-Oscillatory stochastic Galerkin approximation for hyperbolic conservation laws.  Journal of Computational Physics 419, 109663.  doi:https://doi.org/10.1016/j.jcp.2020.109663.
  * Shu [2003] Shu, C.W., 2003.  High-order Finite Difference and Finite Volume WENO Schemes and Discontinuous Galerkin Methods for CFD.  International Journal of Computational Fluid Dynamics 17, 107–118.  doi:https://doi.org/10.1080/1061856031000104851.
  * Shu et al. [2017] Shu, R., Hu, J., Jin, S., 2017\.  A Stochastic Galerkin Method for the Boltzmann Equation with Multi-Dimensional Random Inputs Using Sparse Wavelet Bases.  Numerical Mathematics: Theory, Methods and Applications 10, 465–488.  doi:https://doi.org/10.4208/nmtma.2017.s12.
  * Sullivan [2015] Sullivan, T.J., 2015.  Introduction to Uncertainty Quantification. volume 63.  Springer.  doi:https://doi.org/10.1007/978-3-319-23395-6.
  * Tryoen et al. [2010] Tryoen, J., Le Maître, O., Ndjinga, M., Ern, A., 2010\.  Intrusive Galerkin methods with upwinding for uncertain nonlinear hyperbolic systems.  Journal of Computational Physics 229, 6485–6511.  doi:https://doi.org/10.1016/j.jcp.2010.05.007.
  * Wan and Karniadakis [2005] Wan, X., Karniadakis, G.E., 2005\.  An adaptive multi-element generalized polynomial chaos method for stochastic differential equations.  Journal of Computational Physics 209, 617–642.  doi:https://doi.org/10.1016/j.jcp.2005.03.023.
  * Wiener [1938] Wiener, N., 1938.  The Homogeneous Chaos.  American Journal of Mathematics 60, 897–936.  doi:https://doi.org/10.2307/2371268.
  * Wu et al. [2017] Wu, K., Tang, H., Xiu, D., 2017.  A stochastic Galerkin method for first-order quasilinear hyperbolic systems with uncertainty.  Journal of Computational Physics 345, 224–244.  doi:https://doi.org/10.1016/j.jcp.2017.05.027.
  * Xing et al. [2010] Xing, Y., Zhang, X., Shu, C.W., 2010.  Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations.  Advances in Water Resources 33, 1476–1493.  doi:https://doi.org/10.1016/j.advwatres.2010.08.005.
  * Xiu [2010] Xiu, D., 2010.  Numerical Methods for Stochastic Computations: A Spectral Method Approach.  Princeton university press.
  * Xiu and Hesthaven [2005] Xiu, D., Hesthaven, J.S., 2005\.  High-Order Collocation Methods for Differential Equations with Random Inputs.  SIAM Journal on Scientific Computing 27, 1118–1139.  doi:https://doi.org/10.1137/040615201.
  * Xiu and Karniadakis [2002] Xiu, D., Karniadakis, G.E., 2002\.  The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations.  SIAM journal on scientific computing 24, 619–644.  doi:https://doi.org/10.1137/S1064827501387826.
  * Xiu and Shen [2009] Xiu, D., Shen, J., 2009.  Efficient stochastic Galerkin methods for random diffusion equations.  Journal of Computational Physics 228, 266–281.  doi:https://doi.org/10.1016/j.jcp.2008.09.008.
  * Zala et al. [2020] Zala, V., Kirby, M., Narayan, A., 2020.  Structure-Preserving Function Approximation via Convex Optimization.  SIAM Journal on Scientific Computing 42, A3006–A3029.  doi:10.1137/19M130128X.