Energy Stable and Structure-Preserving Schemes for the Stochastic Galerkin Shallow Water Equations

The shallow water flow model is widely used to describe water flows in rivers, lakes, and coastal areas. Accounting for uncertainty in the corresponding transport-dominated nonlinear PDE models presents theoretical and numerical challenges that motivate the central advances of this paper. Starting with a spatially one-dimensional hyperbolicity-preserving, positivity-preserving stochastic Galerkin formulation of the parametric/uncertain shallow water equations, we derive an entropy-entropy flux pair for the system. We exploit this entropy-entropy flux pair to construct structure-preserving second-order energy conservative, and first- and second-order energy stable finite volume schemes for the stochastic Galerkin shallow water system. The performance of the methods is illustrated on several numerical experiments.


Introduction
The one-dimensional Saint-Venant system of shallow water equations (SWE) is a popular model of water flows where vertical length scales are much smaller than horizontal ones [1].This system in conservative form is given by, where  =  (, ) is the vector of conservative variables; ℎ(, ) is the water height (a mass-like variable) and (, ) is the water discharge (a momentum-like variable).The flux  and source term  are given by, where () is the (assumed known) bottom topography function and  > 0 is the gravitational constant.The system (1.1) is supplemented with initial and boundary data that we omit for the time being.The one-dimensional SWE model (1.1) is a hyperbolic system of partial differential equations (PDE) if ℎ > 0, and hence with the non-zero source , then (1.1) is a nonlinear hyperbolic balance law.Because of this, it inherits the standard challenges in developing numerical methods for such models: solutions generically develop discontinuities in finite time even with smooth initial data, non-uniqueness of weak solutions should be rectified by an implicit or explicit numerical imposition of entropy conditions, and implicit time-integration solvers are challenging to implement due to the nonlinearity [10,27,28].In addition to all this, the SWE has challenges that are somewhat specific to its particular form: positivity of the water height ℎ should be maintained, and numerical schemes should accurately capture near-equilibrium dynamics, which is typically achieved by imposing the well-balanced property [3], i.e., that the PDE equilibrium states are exactly captured at the discrete level.
A more nebulous and hence more frustrating challenge is that of uncertainty in the model.For example, one may have incomplete, partial information about the initial data or the bottom topography function .
In such cases, one frequently models this data as a random variable or process, and hence the solution  to (1.1) is random.We consider the somewhat more simple situation when the input uncertainty is encoded with a finite-dimensional random variable, in which case (1.1) becomes a parametric model (with the input random variables serving as the parameters).Even with this simplification, the parametric or stochastic nature of the solution exacerbates many of the previously described numerical challenges.A particularly successful approach for handling such problems that we will employ is the polynomial Chaos (PC) method, wherein  is approximated as a polynomial function of the input parameters [34,35,44].
The class of non-intrusive PC strategies construct the polynomial by collecting an ensemble of solutions to (1.1) at a collection of fixed values of the parameters.This approach is attractive since it can exploit existing and trusted legacy solvers for (1.1), for which there are several effective choices [4, 9, 14, 23-25, 29, 32, 39-43, 45, 47].However, this approach suffers from the disadvantage that making concrete statements about the quality or properties of the resulting polynomial approximation can be challenging.For example, one cannot guarantee that entropy conditions are satisfied if the polynomial approximation is evaluated away from the parameter ensemble used to construct the approximation.
This paper is concerned with an alternative intrusive approach, the stochastic Galerkin (SG) method for PC approximation, which addresses the parametric dependence in a Galerkin fashion, e.g., by enforcing that certain probabilistic moments of (1.1) vanish.This approach has the potential to provide pathways to mathematical rigor of numerical methods through weak enforcement of the parametric dependence.SG methods transform a parametric model (1.1) into a new non-parametric model of larger system size.Since the new SG formulation is non-parametric, one can apply typical deterministic numerical methods for systems of PDEs to solve the SG problem.Such approaches have shown particular success for modeling parametric dependence in elliptic problems; see, e.g., [8].However, the notable drawback of SG methods when applied to (nonlinear) hyperbolic PDEs is that the new non-parametric SG system need not be a hyperbolic PDE itself, which changes the essential character of the SG system relative to the original system.Despite this challenge, recent work has investigated numerical methods for stochastic Galerkin methods for various types of hyperbolic conservation laws.Some advances include SG-type analysis and algorithms for scalar conservation laws [46] including well-balanced methods [22], Haar wavelet-based SG approaches [20], hyperbolicity preservation through a non-equivalent Roe variables formulation [19], filtering strategies for SG systems [26], limiter-type methods to maintain hyperbolicity [33], hyperbolicity formulations for linear problems [31] or using linearization techniques [38], splitting type approaches for SWE models [7], non-conservative formulations of SG SWE systems [5] and approximate representations through entropic variables [30].
Our focus starts from recent work that has developed an SG formulation for the SWE in conservative form that involves a special SG treatment for the nonlinear, non-polynomial terms [11].This is a particular distinction of our approach: We require no non-conservative formulation, transformation, numerical filtering/limiting, or linearization.Such an approach can be used to develop a well-balanced, hyperbolicity-preserving, and positivity-preserving finite volume method to solve the SG SWE system.The approach forming our starting Table 1.Notation and terminology used throughout this article.

𝐾
Number of terms in a PC expansion ︀  = ( ̂︀ ℎ ⊤ , ̂︀  ⊤ ) ⊤ ∈ R  × R  PC state vector and its components for SG SWE  = (ℎ ⊤ ,  ⊤ ) ⊤ ∈ R  × R  Quantities related to cell averages/reconstructed values, etc. of the PC state vector and its components for SG SWE (, ) Entropy pair for SG SWE ︀  PC vector for velocity  Quantity related to cell averages/reconstructed values, etc. of the PC vector for velocity ︀  , The PC vector for the entropic variable and quantities related to cell averages/reconstructed values, etc. of the entropic variable (•) The operator that maps a PC vector to a PC triple-product matrix ℎ + Energy-conservative flux, 1st-order energy-stable flux, and 2nd-order energy stable flux, respectively  ± , ̃︀  ± Scaled variables point has also been extended to two-dimensional SWE systems [12], but we focus on the single spatial dimension case.

Contributions of this article
We make the following contributions that build on [11]: -We derive an entropy-entropy flux pair for the spatially one-dimensional hyperbolicity-preserving, positivitypreserving SG SWE system derived in [11], see Theorem 3.1.Entropy-entropy flux pairs are the theoretical starting point for proposing entropy admissibility criteria to resolve non-uniqueness of weak solutions.-Using the entropy-entropy flux pair, we devise second-order energy conservative, and first-and second-order energy stable finite volume schemes for the SG SWE, all of which are also well-balanced.See Theorems 4.1-4.3,with the procedure in Algorithm 1.The designed energy conservative and energy stable schemes are stochastic extensions of the schemes developed in [17,18].-We provide numerical experiments that explore the simulation capabilities of the new schemes.To the best of our knowledge, these are the first schemes for any SG SWE system that boast energy stability, the well-balanced property, while also being positivity-and hyperbolicity-preserving.
An outline of this paper is as follows: Section 2 introduces our notation, along with background on PC methods and the SG SWE system from [11].Section 3 provides our entropy-entropy pair construction for the SG SWE system.Section 4 provides the statement of the energy conservative and energy stable schemes that we develop, along with proofs of their theoretical properties, as well as their algorithmic details.Section 5 compiles numerical examples that demonstrate the performance of our scheme.Section 6 gives brief summary of the main results and some future research directions.We summarize our notation in this article in Table 1.
D. DAI ET AL.
In work on the SWE system (1.1) it is common to introduce the water velocity (equilibrium) variable and we also make use of this variable in what follows.

Polynomial chaos expansion
In this section, we briefly review the results and notation for polynomial chaos expansion.More comprehensive results can be found in [13,35,44], etc.
Let  ∈ R  be a random variable associated with Lebesgue density function .Define the function space Assuming finite polynomial moments of all orders for , there exists an orthonormal basis for all , ℓ ∈ N, where  ,ℓ is the Kronecker delta.PC seeks a representation of a random field (•, •, ) ∈  2  in terms of a series of orthonormal polynomials for , (, , ) where ,  are the deterministic spatial and temporal variables, and ̂︀   (, ) are deterministic Fourier-like coefficients.The equation (2.3) holds true for all (, ; •) ∈  2  under mild conditions [15].In practice, a finite truncation of (2.3) is usually considered.Let  be a -dimensional polynomial subspace of  2  , i.e., we let   be an orthonormal basis for  .We make the common assumption that 1 ∈  , and for convenience we assume that, A popular choice for  is the total degree space, but several other options are possible.One choice of -term PC approximation of a random field  in  is the projection of (2.3) onto  : Using the orthogonality of the basis function, the statistics of Π  [] can be expressed in terms of the expansion coefficients.For example, the mean and the variance of Π  [] are given by: A useful lemma is given as follows.
Lemma 2.1.For any two vectors ̂︀ , ̂︀  ∈ R  , The proof is straightforward using (2.7) and (2.8) along with the symmetry of (•).This result is a "commutative" property of the operator (•).For example: For any , ,  ∈ R  , A stochastic Galerkin (SG) formulation of a -parameterized PDE corresponds to making the ansatz that the state variable lies in the space  , and projecting the PDE residual onto the same space.Straightforward applications of this procedure to (nonlinear) hyperbolic PDEs typically do not result in hyperbolic SG formulations.

Hyperbolic-preserving stochastic Galerkin formulation for shallow water equation
In [11], we have derived a hyperbolicity-preserving stochastic Galerkin formulation for the shallow water equations.We briefly recall the results in this section.
We make the ansatz that ℎ,  lie in the polynomial space  , and use these to formulate a -variable Galerkin projection of the SWE.We make a special choice of how the D. DAI ET AL.
Galerkin projection of the nonlinear, non-polynomial term () 2 /ℎ is truncated, which results in a new (stochastic Galerkin) system of balance laws whose state variables are the expansion coefficients in (2.11) [11]: Here, ̂︀  := ( ̂︀ ℎ ⊤ , ̂︀  ⊤ ) ⊤ ∈ R 2 , where ̂︀ ℎ, ̂︀  are each length- vectors of the expansion coefficients in (2.11).The flux and the source terms are, cf. (1.2).The flux Jacobian, written in  ×  blocks, is given by We have introduced the term which we view as the vector of the PC coefficients of the -velocity  introduced in (2.1), and is well-defined if ( ̂︀ ℎ) is invertible.
The deterministic SWE are hyperbolic if the water height ℎ > 0; there is a natural extension of this property to the SG SWE.
Theorem 2.1 ( [11], Thm.3.1).If the matrix ( ̂︀ ℎ) is strictly positive definite for every point (, ) in the computational spatial-temporal domain, then the SG formulation (2.12) is hyperbolic.This is proven by identifying a stochastic extension of the known eigenvector matrix for the deterministic SWE flux Jacobian   , and using this to show that  ̂︀   ̂︀  is similar to a symmetric matrix and hence (2.12) is hyperbolic [11].

An entropy-entropy flux pair for SG SWE systems
The formulation (2.12) will be considered in what follows.Our goal will be to derive entropy-entropy flux pairs for these formulations.The first step is for us to recall a known entropy-entropy flux pair for the deterministic SWE system.

Entropy-entropy flux pairs for deterministic shallow water equations
It is well-known that solutions to systems of conservation/balance laws can develop shock discontinuities in finite time for generic initial data.Therefore, weak solutions, i.e., solutions in the sense of distributions, are usually considered.However, weak solutions are not necessarily unique, and to mitigate this issue an additional entropy admissibility criteria is imposed [2,10] to identify the physically meaningful solution.
For a general balance law in one space dimension its entropy-entropy flux pair (( ), ( )) satisfies a companion balance law where the entropy ( ) is a scalar function that is convex in  , and  is an entropy flux function.In order to be consistent with the original balance law for smooth  , the entropy-entropy flux pair (, ) should satisfy the following compatibility condition, which is simply the condition ensuring that multiplying (3.1) by   recovers (3.2) when solutions are smooth.In the case of  ≡ 0 and (, ) = (( ), ( )), equation (3.3) is the usual entropy condition for conservation laws.For a general system of balance laws in several spatial dimensions, an entropy-entropy flux pair need not exist.However, for a hyperbolic system of balance laws emerging from continuum physics, the companion balance law (3.2) is usually related to the Second Law of thermodynamics, and the total energy of the system often serves as the entropy function.A variety of examples can be found in Section 3.3 from [10].For the deterministic SWE system in (1.1), the total energy [17] is where we recall that  is the velocity defined in (2.1).For any smooth solution  , a direct calculation yields, where This, along with the fact that   is convex in  , establishes that (  ,   ) is a valid entropy-entropy flux pair for (1.1).For (weak) solutions with shocks, the entropy admissibility criteria is that energy should dissipate in accordance with a vanishing viscosity principle, In what follows we will identify entropy-entropy flux pairs for the SG SWE model.This amounts to verifying that (i) such a pair satisfies the companion balance law (an equality for smooth solutions) and (ii) that the entropy function is convex in the state variable.

An entropy-entropy flux pair for the one-dimensional SG SWE
This section is dedicated to identifying an entropy-entropy flux pair for the SG system (2.12).In this section, we will return to the notation ̂︀  (containing PC expansion coefficients) for the derivation of an entropy-entropy flux pair for the SG system.Our main result in this section is the following entropy entropy-flux pair for the one-dimensional SG SWE: Theorem 3.1.Define the function, and also the flux function, If ( ̂︀ ℎ) > 0, then (, ) is an entropy-entropy flux pair for the one-dimensional SG SWE (2.12).
Recall that ̂︀  above is defined in (2.15), and contains PC expansion coefficients for the velocity  defined in (2.1).In the absence of uncertainty, equation (3.8a) reduces to the deterministic total energy (3.4).The rest of this section is devoted to proving Theorem 3.1, which amounts to showing that, if ( ̂︀ ℎ) > 0, then  is convex in ̂︀  and (, ) satisfy the companion balance law, for smooth solutions ̂︀  .Note that for non-smooth solutions, equation (3.9) holds with = replaced by ≤.We prove Theorem 3.1 with three intermediate results.Our first result is a technical condition that facilitates later computations.
Proof.If () is a -parameterized matrix, then for any  at which  is invertible, Applying this to , we have, and hence, proving the desired relation for Proof.Using the definition (2.15) of ̂︀ , note that, and therefore in particular, We will show that this Hessian is positive definite.Clearly we have, Using the previous lemma, we can directly compute, which in turn implies, and finally, Hence, the Hessian of  1 is, A direct computation of the quadratic form associated to this Hessian using an arbitrary vector Finally, combining the above with (3.15) and (3.16) yields, which is non-negative since ( ̂︀ ℎ) is positive-definite.Therefore,  is convex, as desired.In addition, since the above expression vanishes if and only if  1 =  2 = 0, then  is also strictly convex.
The final piece needed to prove Theorem 3.1 is to establish that the entropy function  along with the flux function  defined in (3.8b) satisfy the companion balance law.

Lemma 3.3 ((𝐸, 𝐻) satisfy the companion balance law). When ̂︀
is a smooth function, the pair (, ) defined in (3.8) satisfies Proof.The compatibility condition we seek to show, equivalent to (3.18), is, cf. (3.3).To proceed we split both entropy functions into two pieces: From (3.14), (3.16), and (3.17), we have already computed the gradient of : Combining these expressions with the flux Jacobian in (2.14) and the source term in (2.13) yields, Note then that if we are able to show, then the expressions (3.22) are equivalent to (3.19).Therefore, we are left only to show (3.22b).A direct computation with (3.21) and (2.14) yields, On the other hand, we have the expressions, and using these to compute 1  ̂︀  shows that (3.22b) is true, completing the proof.
The proof of Theorem 3.1 is complete: Lemmas 3.2 and 3.3 imply that (, ) as defined in (3.8) are an entropy-entropy flux pair for (2.12).
Remark 3.1.The quantities, are called the entropy variable and stochastic energy potential, respectively.These variables serve important roles in the construction of the energy conservative and the energy stable schemes that we develop later.

Well-balanced energy conservative and energy stable schemes for the SG SWE
In this section, we present several well-balanced energy conservative and energy stable numerical scheme for the SG SWE.The schemes designed below are stochastic extensions of the schemes developed in [17].Our entropy-entropy flux pairs developed in Section 3 will be crucial ingredients for energy conservative and energy stable schemes for the SG formulation (2.12)- (2.15).
We also need to specify the well-balanced property we are interested in: By "well-balanced", we mean that the scheme can preserve the stochastic "lake-at-rest" state exactly at the discrete level.Definition 4.1 (Well-Balanced SG SWE Property, [11]).We say that a solution (ℎ  ,   ) to (2.12) is wellbalanced if it satisfies the stochastic "lake-at-rest" solution, where () is a random scalar depending only on , Π  corresponds to a polynomial truncation, cf.(2.5), and subscripts  refer to the (stochastic) discrete solutions in the subspace  .In terms of our previous notation for  -expansion coefficients, equation (4.1) is equivalent to the following vector equation where  is the spatial domain and  is the terminal time.
We emphasize that even without introducing the lake-at-rest definition (4.1), the vector equation (4.2) itself is a steady state of the SG system (2.12).

Energy conservative schemes
We consider the semi-discrete form for FV schemes for (2.12) over a uniform mesh in the  variable: d is a discretization of the source term, which we will design below to be well-balanced.To reiterate our notation: normal typeset capital letters (sometimes with "hat" notation) refers to degrees of freedom associated to discretizing only the stochastic variable , i.e., ( ̂︀  , ̂︀ ℎ, ̂︀ , ̂︀ ).Boldface notation with subscripts  refers to degrees of freedom associated to a subsequent discretization of the spatial variable  over cell ℐ  , i.e., (  , ℎ  ,   ,   ).We define the discrete velocity variable   in a manner analogous to (2.15): Discrete entropic quantities are derived from the discrete conservative variables   and velocity variable   .I.e., the following are direct generalizations of the definition of ( ̂︀  ) in (3.8a), and of ( ̂︀  , Ψ) in (3.23): ) We now introduce some notation that is used in [17] for averages and jumps at cell interfaces: where   is the cell average over ℐ  .The expressions above are equivalent to, and all these expressions are valid regardless of the size of  (e.g., both row and column vectors are allowed).
We will require some additional technical results for interfacial averages and jumps.
We now make particular definitions for energy conservative and energy stable schemes for one-dimensional systems of balance laws.To provide context, with no source terms (i.e.,   = 0) then spatial discretizations of the form (4.3) are called conservative schemes since they imply, and in particular with periodic boundary conditions, then this implies that the cumulative amount of ̂︀  in the system is constant in time 1 .
To translate this concept to the notion of an energy conservative scheme, note that an entropy-entropy flux pair (, ) introduced in Section 3 is explicitly a function of the state ̂︀  and inputs in the source term (here, ︀ ).Hence, the semi-discrete form (4.3) can be transformed into a semi-discrete form for the companion balance law (3.9).Then we call (4.3) energy conservative if it implies a conservative scheme for the companion balance law that describes the evolution of the entropy (energy).
Definition 4.2 (Energy conservative and energy stable schemes).Suppose that the system of balance laws (2.12) has an entropy-entropy flux pair (, ) where ( ̂︀  ) can be interpreted as energy for the system.Then the semi-discrete FV scheme (4.3) is an Energy Conservative (EC) scheme if it can be rewritten as the following semi-discrete form for the evolution of the numerical cell averages   of : where ℋ +1/2 is some numerical entropy flux at the interface location  =  +1/2 .The scheme (4.3) is called an Energy Stable (ES) scheme if Note that the definitions above are cell-wise conditions that are stronger than a global condition such as (4.9).

An EC scheme for the SG SWE
In this section we present an EC scheme for the one-dimensional SG SWE system (2.12).We use the conservative scheme (4.3), with the following choices of flux and source terms: ) Above, the interfacial averages  +1/2 are computed as defined in (4.6).Our main result for this scheme is as follows.
Theorem 4.1 (EC Scheme).Suppose the bottom topography function  is independent of time.Consider the semi-discrete scheme (4.3) for the SG SWE system (2.12).Suppose that the flux and source terms are selected as in (4.12).Then, this is a well-balanced EC scheme with local truncation error (∆ 2 ).
The remainder of this section is devoted to the proof, which requires some intermediate steps.First, we show that   is a well-balanced choice for the source term discretization.Lemma 4.2.Suppose   is chosen as in (4.12b).If the bottom topography  is independent of time, then (4.3) is a well-balanced scheme in the sense of Definition 4.1.

Proof. Given initial data
the well-balanced property with time-independent bottom topography (see Def. 4.1) requires that, for every , We first notice that, Note that with initialization (4.13), then   = 0, and hence  +1/2 = 0. Therefore the semi-discrete scheme (4.3) with the flux and source terms in (4.12) yields, Note that the implicit constants hidden in the asymptotic notation above depend on the maximum singular value of ( ̂︀ ℎ  ()) and the minimum singular value of ( ̂︀ ℎ( +1/2 )).
The final result we need is a sufficient condition for a numerical flux to result in an EC scheme.
We now have all the ingredients necessary to prove Theorem 4.1.
Proof of Theorem 4.1.Lemmas 4.2 and 4.3 verify that the scheme is well-balanced and second-order.We therefore need only show that it is EC.To do this, we must verify the condition in Lemma 4.4.We accomplish this with direct computation: which verifies (4.17), and hence Lemma 4.4 is applicable, showing that this is an EC scheme.

A first-order ES scheme
The scheme determined by (4.12) numerically preserves the energy of the PDE system (1.1).However, it may lead to spurious oscillations since the energy should dissipate in the presence of shocks.The issue can be resolved by introducing appropriate numerical viscosity [16-18, 36, 37].Our numerical diffusion operators are a straightforward stochastic extension of the energy-stable diffusion operators proposed in [16,17].
For context of the approach, the introduction of a traditional Roe-type diffusion for a conservation law involves augmenting an EC flux as follows: where  Roe is a positive semi-definite matrix defined through a diagonalization of the interfacial flux Jacobian at a Roe-averaged state: Then the semi-discrete scheme (4.3) using the numerical flux ℱ +1/2 = ℱ RD +1/2 would behave like, where  is a positive-definite matrix and ̂︀   is the second spatial derivative of the PDE state variables introduced in equation (2.12), and hence this introduces diffusion into an EC scheme.While the above approach works in terms of adding a diffusion-like term, a convenient way to ensure energy stability is to employ a numerical diffusion term that operates on the entropic variables  instead of the conservative variables  : where  ES is a positive definite matrix that will be identified in a Roe-type way from the two adjacent states   and  +1 at the cell interface  =  + 1  2 .The term   is as given in (4.5b), and is a second-order approximation to the cell-average of the entropy variable ̂︀  .We are interested in the Roe-type energy-stable operator defined as, where the matrices  and Λ are matrices from the eigendecomposition of the flux Jacobian (2.14) evaluated at a Roe-type average state: Note in particular that  +1/2 ̸ = (ℎ +1/2 ) +1/2 , so that ̃︀  +1/2 ̸ =  +1/2 .The focal scheme of this section uses the numerical flux (4.21),where  is given by the Roe-type diffusion matrix introduced above, where we refer to this scheme as "ES1" because we will show it is first-order accurate.Our main result for this scheme is as follows.
Theorem 4.2 (ES1 scheme).Consider the finite volume scheme (4.3) with source term (4.12b) and diffusive numerical flux (4.21), selecting the diffusion matrix as, The resulting scheme is a first-order, well-balanced ES scheme.
Proof.We omit some details that are similar to the proof of Theorem 4.1.We have already established in Theorem 4.1 that ℱ EC +1/2 is second-order accurate.That this ES1 scheme is first-order is direct from the definition of   in (4.5b), resulting in the approximation which implies that the diffusive augmentation in (4.21) commits a first-order local truncation error.To establish that this scheme is well-balanced, we assume the stochastic lake-at-rest initial data (4.13), and this coupled with the definition of   in (4.5b) implies  +1/2 = 0. Since the EC flux and source are well-balanced (Lem.4.2), then this implies that this ES1 scheme is also well-balanced.
Finally, we seek to show the ES property.We define the ES1 energy flux, with ℋ +1/2 as defined in (4.18).As in Lemma 4.4, we multiply (4.3) by  ⊤  ; after manipulations that are similar to those in the proof of Lemma 4.4, we have, )︁ .
Since  ES1 +1/2 is positive semi-definite, then this scheme satisfies (4.11), and hence is an ES scheme.

ES1 diffusion vs. Roe diffusion
We provide in this section a result that motivates and justifies our particular form of the ES1 diffusion modification defined in (4.24) and (4.25).This result states that if the bottom topography function vanishes (i.e., we are in the specialized case of a conservation law), then our chosen Roe-type ES1 diffusion in (4.22) and (4.23) coincides with a standard Roe-type diffusion term.Hence, in specialized scenarios our diffusive augmentations using entropic variables are equivalent to more standard Roe-type diffusion.Proposition 4.1.Define the Roe diffusion matrix as in (4.20), but using the flux Jacobian evaluated at ̃︀ where we have evaluated the flux jacobian at ̃︀  +1/2 instead of at  +1/2 .Assume   = 0 for all  ∈ [ ].Then, Proving this result requires some setup: Under the assumptions of Proposition 4.1 we consider the SG SWE (2.12) with flat bottom, i.e., ̂︀  = 0, together with entropy  flat ( ̂︀  ) = 1 2 (̂︀ ) ⊤ ̂︀ +  2 ‖ ̂︀ ℎ‖ 2 and the entropy variables,

.28)
Our main tool will be some results of the proof of Theorem 3.1 in [11]; in particular, while we have provided the flux Jacobian for this system in (2.14), we will need the explicit similarity transform that accomplishes its symmetrization.
where  is the symmetric matrix, and The second lemma reveals the relation between the cell interface jump of  flat (the spatial approximation corresponding to the cell-averaged entropy variable ̂︀  flat in (4.28)) and  .
Lemma 4.6.Recall the definition of ̃︀  + 1 2 in (4.23): which is an intermediate state defined by the arithmetic average of ℎ and  across the cell interface  =  + 1 2 .Denote,  flat to be the corresponding spatial approximation of the cell-averaged entropy variable defined in (4.28).Then, where )︃ .
Proof of Proposition 4.1.Let  + 1 2 be the symmetric matrix defined in (4.29) evaluated at  + 1 2 , and be its eigenvalue decomposition.Then, is an eigendecomposition of the Jacobian matrix  ̂︀ , where we have used the fact that  −1 due to the symmetry of  + 1 2 .The Roe-diffusion operator evaluated at the location ̃︀  + 1 2 as indicated in (4.26) is then given by, Therefore, (4.36)

A second-order ES scheme
To develop a second-order accurate energy-stable scheme, we use jump operators with (∆ 2 ) accuracy.A natural choice is to use the jumps obtained by non-oscillatory second-order reconstructions of the entropy variable.However, attaining a provable energy-stable scheme requires the more subtle reconstruction procedure in [18] that we follow.The new idea for second-order diffusion is to use reconstructions in order to compute jumps.To that end, we let  +  and  − +1 be second-order reconstructions from the right and left, respectively, of the entropy variable  () at location  =  +1/2 .We will describe later in this section how these reconstructions are computed.
Assuming we have these reconstructions in hand, we can compute second-order accurate jumps of the entropy variables: The overall scheme is similar as the previous section, but uses a second-order diffusive augmentation of a conservative flux, D. DAI ET AL.
We choose the matrix  ES2 as for the ES1 scheme, where we recall that the eigendecomposition matrices  , Λ are computed from the Roe-type average of the flux Jacobian, cf.(4.22), (4.23).One could alternatively select  ES2 by using second-order reconstructions of  as input to , e.g., for some second-order reconstructions  ±  .What remain is to describe how  ±  are computed in a way that ensures the energy stable property.The main idea is to design  ±  through a second-order reconstruction of scaled (transformed) versions of the entropy variables: where the matrices  ±1/2 are as in (4.39).Once these have been computed, we perform a second-order total variation-diminishing (TVD) reconstruction on the  variable at the interfaces: where ∘ is the Hadamard (elementwise) product on vectors, and  ±  are difference quotients, where ⊘ is the Hadamard (elementwise) division between vectors.We select the function  to be the minmod limiter, which operates elementwise on vector inputs.Note that other slope limiter functions  may be selected, but minmod is the only valid limiter in this context that also satisfies the TVD property ( [18], Sect.3.4).Finally, the desired reconstructions for  ±  are defined by inverting the -to- map, The full scheme has now been described, and satisfies the following properties.We focus the remaining discussion in this section on sketching the proof of the above result.The second-order property results from the fact that the jumps are computed using second-order accurate reconstructions; the well-balanced property can be proven in exactly the same way as is done for the ES1 scheme in the proof of Theorem 4.2.To show the ES property, we exercise one of the major results in [18] that we reproduce below.Lemma 4.7 ([18], Lem.3.2).For each , if there exists a positive diagonal matrix Π +1/2 ≥ 0 such that the second-order jump satisfies, then the scheme (4.3) with flux term ℱ +1/2 = ℱ ES2 +1/2 is an ES scheme.
Hence, showing the ES property for our scheme only requires us to establish (4.44).To accomplish this, note that the definition (4.41) implies, I.e., we have, and in particular Π +1/2 is a diagonal matrix and positive semi-definite since 0 ≤ () ≤ 1.Since the jump operators ⟨⟨•⟩⟩ and • are linear in their arguments, then combining (4.45) with the relations (4.40) and (4.43) that connect   and ̃︀   to   and  ±  yields the relation (4.44) with a positive-definite diagonal matrix Π +1/2 .Hence, this is an ES scheme, and completes the proof of Theorem 4.3.
Finally, we remark that the implementation of the diffusion term in the ES2 flux (4.38) does not require explicit construction of  ±  .I.e., we have, and hence one need only compute ̃︀  ±  in order to directly evaluate the diffusion part of the ES2 flux.

Algorithmic details
Our overall scheme is the semi-discrete form (4.3), which we pair with a numerical time-stepping scheme.We provide pseudocode in this section that describes a fully discrete SG SWE time-stepping algorithm.This full pseudocode introduces some additional details for the scheme that were devised in [11], many of which are based on standard procedures used in schemes for deterministic SWE models [23].We very briefly describe these additional details in the coming sections; more comprehensive discussion can be found in [11].The full algorithmic pseudocode is given in Algorithm 1.

Velocity desingularization
Computing   requires inversion of the matrix (ℎ  ), which is assumed (and enforced in the scheme) to be symmetric and positive-definite.However, this matrix may be ill-conditioned.To ameliorate numerical artifacts associated with this ill-conditioned operation, we employ a desingularization procedure, introduced for the deterministic SWE in [25].We describe here the stochastic variant of the desingularization procedure, proposed in [11].If (ℎ  ) has the eigenvalue decomposition, where   > 0 are the eigenvalues of (ℎ  ), then the desingularization process approximates (ℎ  ) −1   by regularizing the matrix inverse procedure: where  > 0 is a small constant; we choose it to be  = ∆.Note that if   ≥  1/4 , then ̃︀   =   , and hence regularization is performed only in the presence of small eigenvalues.Compared to (4.4).This procedure to compute   is a stabilized way to compute velocities.
For scheme consistency, if the desingularization above is activated, then we recompute the discharge variable: This update for   may affect global conservation of  through a ( 2 /  ) ∼ (∆ 2 /  ) augmentation; however, this desingularization process is constructed to minimize the numerical effects of such desingularization while simultaneously restoring well-behaved and stable computation.Without such a desingularization procedure, the scheme can quickly become numerically unstable for small water heights.

Hyperbolicity preservation
The SG SWE PDE (2.12) is hyperbolic and has an entropy pair if ( ̂︀ ℎ) > 0, i.e., Theorems 2.1 and 3.1, respectively.To ensure this holds at the discrete level, we require the condition (ℎ  ) > 0 for every cell .To enforce this, we employ ( [11], Thm.3.4, Cor.3.5), which state that a sufficient condition for (ℎ  ) > 0 is that for every  = 1, . . .,  , where {  }  =1 is a nodal set in R  for a positive-weight quadrature rule having sufficient accuracy relative to the -polynomial space  defined in (2.4).The functions   are the basis of  in (2.4) for which ℎ  are coordinates.The function ̂︀ ℎ  () is the SG SWE approximation to the ℐ  -cell average of ̂︀ ℎ(, , ) at the current time.Hence, the computational vehicle we use to enforce hyperbolicity of the underlying PDE in our scheme is to enforce the above positivity-type condition on the ℎ  variable.

Positivity-preservation
We enforce the positivity condition (4.48) by restricting the timestep size.We assume that the current time value of ℎ  satisfies (4.48).If forward Euler with a stepsize ∆ is used to discretize (4.3), then (4.48) is true at the next time step if, where ̂︀ ℱ ℎ +1/2 (•) is the SG approximation of the ̂︀ ℎ-variable flux: Hence, we enforce positivity preservation by ensuring a small enough timestep so that the positivity condition (4.48) is respected globally over all spatial cells.We must also restrict ∆ to satisfy the wave speed CFL condition; see equation (4.16) of [11].

Adaptive time-stepping
The time step restriction (4.49) works for forward Euler time-stepping.To extend this to a higher-order temporal scheme, we employ a third-order strong stability-preserving scheme, which is a convex combination of forward Euler steps [21].However, the intermediate stages of a(ny) time-stepping scheme need not obey the positivity-preserving property, even if ∆ is chosen to obey the condition (4.49) determined at the initial step.Legendre polynomials.Our final example in Section 5.5 employs a two-dimensional vector , and the details of that experiment are provided in that section.
Instead of visualizing the conservative variable ℎ corresponding to water height, we will plot the water surface , defined as  = ℎ + , with the bottom topography  superimposed on the same graph; plots of (, ) are more physically interpretable than directly plotting the water height ℎ.

Accuracy test for the EC scheme on smooth initial data
We start by examining the accuracy of the proposed EC scheme (4.3), (4.12a), (4.12b) for the SG system (2.12) on a test with a smooth initial data and a smooth solution (we select a final time to ensure that the solution is smooth).On a smooth solution, the system satisfies the energy equality (3.5) and the EC scheme is constructed to satisfy the discrete version of the energy equality.We consider a smooth stochastic water surface, ℎ(, 0, ) + (, 0, ) = (, 0, ) = 1.1 + 0.1 −2 + 0.001 −10 sin(cos(2)) , (, 0, ) = 0.1, ( with a flat bottom (, , ) ≡ 0. For all tests in this section we compute the solution to a final time  = 0.0025 over the physical domain  ∈  = [−1, 1].The profile of the water surface and discharge computed by the EC scheme at a final time on a mesh of   = 3200 elements with  = 4 is shown in Figure 1 (water surface (left) and discharge (middle)).The right plot in Figure 1 illustrates that the EC scheme numerically preserves energy on this test -a very small numerical error of (10 −12 ) in the energy is due to time discretization.We next test spatial accuracy of the EC scheme against a reference solution computed by the EC scheme on mesh of size   = 3200 with  = 4, Figure 1.We illustrate convergence of the scheme using the water height ℎ by computing the error between the reference solution and the numerical solution using an  1 norm in physical space and an  2 norm in parameter space, where we recall that  is the physical domain and R  is the stochastic domain, ℎ  is the numerical water height PC solution to the SG SWE equation on a mesh of size   , and ℎ ref is the -term (fixed  = 4) PC reference solution.The results for this spatial convergence test are presented in Table 2a and show second order convergence for the smooth solution, as expected.
Our next test investigates convergence with respect to the number of PC terms  for one-dimensional parameter space,  = 1.For infinitely smooth solutions, we expect spectral (often exponential) convergence.For this test, the reference solution corresponds to   = 6400 with  = 25, and again illustrate convergence of the scheme using water height ℎ by computing the error between the reference solution and the numerical solution Table 2. Results for Section 5.1: EC order of convergence in physical space (A) for the problem (5.2) with smooth initial data.Error is computed using (5.3a).Order of convergence in stochastic space is presented in (B).Error is computed using (5.4) where ℎ  is the water height corresponding to a solution using  PC terms.We observe spectral convergence of EC scheme in Table 2b, faster than any fixed polynomial order, suggesting this solution has very high regularity, and confirming the expected numerical convergence rate in .
From Figure 2, similar to the results presented in Figures 1 and 4 of [17], we observe that the water surface with uncertainty develops a leftward-going rarefaction wave and a rightward-going shock.Similar to [17], EC computes such solutions accurately, but at the expense of large post-shock oscillations as observed on Figure 2 (right plot).These oscillations are expected since the EC scheme preserves energy on such solutions, and hence energy is not dissipated across the shock as it should.We also demonstrate in Figure 3 (middle and right plots) numerical energy conservation for the EC scheme.We note that the energy conservation errors due to time discretization are reduced significantly by decreasing the time step/CFL constant (right figure), similar to the results reported in Figure 1 in [17].The presented results in Figures 2 and 3 (left figure) also illustrate that ES2 produces less smearing than ES1 at both the rarefaction and the shock waves.The schemes ES1 and ES2 are both designed to dissipate energy which is also confirmed by the numerical results in Figure 3 (middle plot), with the energy dissipation in ES2 being lower than with the ES1 scheme.In addition, the numerical results seem to indicate that the ES2 scheme is better able to capture large variance spikes compared to the ES1 scheme.Finally, the employment of the the numerical diffusion operators in ES1 and ES2 schemes removes oscillations present in the numerical solution using the EC scheme.The observed results are also in agreement with deterministic model numerical results reported in [17].

Stochastic bottom topography
Next, we consider the shallow water system with deterministic initial conditions, (, 0) = {︂ 1  < 0 0.5  > 0 , (, 0) = 0,  This test example was presented previously in [11].Initially, the highest possible bottom barely touches the initial water surface at  = 0, see Figures 4-6.We use a uniform grid size   = 400, 800, 1600 over the physical domain  ∈ [−1, 1], and compute up to time  = 0.0995.(Immediately after this time, the EC scheme fails for   = 400 due to spurious oscillations near sharp gradients of the solution.)In Figures 5 and 6 we compare only performance of ES1 and ES2 at  = 0.0995 since EC fails on those meshes even earlier.Again, the numerical results indicate that the ES2 scheme can more easily resolve large, spatially-concentrated variance values compared to the ES1 scheme, but under mesh refinement both schemes converge to similar numerical solutions.In Figure 7, we show the numerical solution obtained using ES1 and ES2 at the final time  = 0.8 and on a mesh of size   = 800.For both schemes, the 99% confidence region of the water surface stays above the 99% confidence region of the bottom function in Figure 7, and both methods produce similar numerical solutions.The presented results are comparable to the results in Section 5.1 of [11].In Figure 8, we again observe   as expected that the EC scheme numerically conserves energy, while ES1 and ES2 dissipate energy with larger dissipation produced by ES1 method.(5.7)

Perturbation to lake at rest
The test is from [7] and is similar to the deterministic tests of the perturbation of lake at rest solution, for example to the one presented in [17].We start by illustrating the accuracy of the ES1 and ES2 schemes on this problem with discontinuous initial data and nonsmooth solution, in which case we expect the system to dissipate energy.We compute the numerical    solution to a final time  = 0.8.We test the spatial accuracy of ES1 and ES2 schemes against reference solutions computed by ES1 and ES2 schemes, respectively, computed on a mesh of size   = 3200 with  = 2.As before, we illustrate convergence of the scheme using the water height ℎ by computing the error between the reference solutions and the numerical solutions using (5.3a).
From Table 3 we observe that both ES1 and ES2 schemes perform well on problems with nonsmooth solutions (the order of convergence is limited by the regularity of the solution rather than the formal order of the scheme).ES1 produces somewhat higher than the guaranteed first order of convergence, while ES2 produces slightly lower than second order of the convergence, which is expected on such nonsmooth problems.The ES2 scheme delivers overall better resolution and accuracy than the ES1 scheme, as can be seen from the table and the plots of the solutions at different times in Figures 9 and 11.
From the results in Figures 9-11, we make conclusions similar to previous sections: both ES1 and ES2 capture small stochastic perturbations of the lake at rest solution quite well (with both leftward-and rightward-going waves present in the numerical solutions).The first order ES1 scheme exhibits much more dissipation in the left-and the right-going waves than the ES2 scheme, which produces a more accurate solution, as shown in Figures 9, 10 (left and middle figures), and 11.The results of the EC scheme is also shown in Figure 10 (right figure).The EC scheme resolves the left-and the right-going water waves with heights higher than in both ES1 and ES2 methods, but again there are oscillations present near both waves in the EC numerical solution since    EC does not dissipates energy across shocks.The relative energy change for this example produced by EC, ES1 and ES2 methods is illustrated in Figure 12.The presented results are also comparable to the results reported in [17] and in [7].(5.9)In this system we consider the two-dimensional random variable  = ( 1 ,  2 ), where  1 and  2 are modeled as independent and identically distributed random variables over [−1, 1] having a Beta distribution with parameters (, ) = (1, 3).The plot of the bottom mean and quantiles as defined by  in (5.9) is presented in Figure 13 (left).To simulate the SG SWE system, we use  1 = 3 and  2 = 5 polynomial terms for parameters  1 and  2 , respectively, corresponding to a tensor-product space of polynomials having dimension  = dim  =  1  2 = 15.We simulate up to time  = 0.8 with   = 400 elements for  ∈ [−1, 1].The numerical results are plotted in Figure 14, showing the water surface and the discharge develop more complex structures than in the previous example (5.6) and (5.7) in Section 5.4, and the relative energy change shown in Figure 13 (right) is of a comparable order of magnitude.In addition, we draw conclusions similar to previous sections: both ES1 and ES2 capture small stochastic perturbations of the lake at rest solution quite well (with both leftward-and rightward-going waves present in the numerical solutions).The first order ES1 scheme exhibits more dissipation in the left-and the right-going waves than the ES2 scheme, which produces a more accurate solution, as shown in Figure 14.

Conclusion
In this work we derived an entropy-entropy flux pair for the spatially one-dimensional hyperbolicitypreserving, positivity-preserving SG SWE system developed in [11].Such entropy-entropy flux pairs are the theoretical starting point for proposing entropy admissibility criteria to resolve non-uniqueness of weak solutions.Using the proposed entropy-entropy flux pair, we designed second-order energy conservative, and first-and second-order energy stable finite volume schemes for the SG SWE.The proposed schemes are also well-balanced.We provided several numerical experiments to illustrate performance of the methods.As part of future research, we plan to extend such methods to models in two spatial dimensions, to explore alternative constructions of diffusion operators, and to investigate other reconstruction approaches for the entropy variables.

1 .
Results for Section 5.1.Left: water surface and middle: discharge.Right: relative energy change in EC over time.Mesh size   = 3200 and PC basis functions  = 4.

Figure 7 .
Figure 7.Comparison of the results for Section 5.3 using different schemes.Top: water surface.Bottom: discharge.Left: ES1, and right: ES2.Mesh   = 800 with  = 9 at the final time  = 0.8.

Table 3 .
Results for Section 5.4: ES1 and ES2 order of the convergence in space for nonsmooth solution (5.6) and (5.7).Error is computed using (5.3a).