Entropy stability and well-balancedness of space-time DG for the shallow water equations with bottom topography

We describe a shock-capturing streamline diffusion space-time discontinuous Galerkin (DG) method to discretize the shallow water equations with variable bottom topography. This method, based on the entropy variables as degrees of freedom, is shown to be energy stable as well as well-balanced with respect to the lake at rest steady state. We present numerical experiments illustrating the numerical method.


Introduction.
The shallow water equations model many phenomena of interest to meteorology and oceanography.In one space-dimension, the shallow water equations with variable bottom topography are given by, h t + (hu) x = 0, (hu) t + 1 2 gh 2 + hu 2 x = −ghb x . ( Here h is the water height, u is the depth-averaged velocity, g the gravitational constant, and b is the bottom topography.The two-dimensional version of the equations is provided in the appendix.The shallow water equations with bottom topography are an example of nonlinear systems of balance laws of the form, by defining Here and in the following, we drop the index k in the one-dimensional case for simplicity. It is well-known that solutions of balance laws (or conservation laws, when the source term S(U) ≡ 0) can be discontinuous due to the appearance of shock waves (hydraulic jumps for shallow water equations), even when the initial data are smooth.Hence, the solutions to such equations are sought in the weak sense [3].However, weak solutions are not unique and one needs to impose additional admissibility criteria, in the form of entropy conditions, to recover uniqueness [3].In general, entropy conditions are of the form, with the above entropy inequality being satisfied in the sense of distributions.
In the specific case of shallow water equations with variable bottom topography, it is well-known that the total energy, S(U) = 1 2 hu 2 + 1 2 gh 2 + ghb, Q(U) = 1 2 hu 3 + gh 2 u + ghbu, serves as the entropy with the energy flux Q.
Furthermore, the entropy variables V = S U (U), with the specific form, play an important role in the entropy stability analysis.In particular, we can rewrite (2) in entropy variables as using the notation F k (V) = F k (U(V)) and S(V) = S(U(V)) for simplicity.
1.1.Numerical methods.The design of efficient numerical methods for approximating systems of conservation (balance) laws is fairly mature.Among the popular discretization frameworks are the finite volume (conservative finite difference) methods, based on (approximate) Riemann solvers [12].High-order of spatial accuracy can be attained by employing suitable non-oscillatory piecewise polynomial reconstructions such as TVD, ENO and WENO.An attractive alternative is provided by the discontinuous Galerkin method.
In the specific case of shallow water equations, a further challenge is provided by the fact that most situations of interest are modeled as small perturbations of a steady state, the so-called lake (ocean) at rest state [13] given by, Therefore, one needs to design numerical methods that preserve this steady state.If a numerical method only approximates this steady state to truncation error, then it will not be able to resolve the very small amplitude waves (perturbations of steady state) that are of interest, except on very fine meshes.Hence, the goal has been to design well-balanced schemes, i.e. numerical approximations of the shallow water equations that preserve a discrete version of the lake at rest steady state [13].
A large number of well-balanced schemes have been designed for the shallowwater equations over the past decade or so.A (very incomplete) list of references is [1,9,10,5,2,11,14] and references therein.The basic idea behind most of these papers is to modify the numerical fluxes by a hydrostatic reconstruction and introduce a source discretization to balance the flux difference.
Unfortunately, very few rigorous stability results exist in the context of wellbalanced schemes for the shallow water equations with variable bottom topography.
Given the fact that energy stability appears to be the most natural stability framework, it is natural to seek a well-balanced scheme that is also entropy stable.To the best of our knowledge, this issue is first tackled in [4], where the authors designed an energy stable scheme for the shallow water equations with bottom topography.Additionally, this scheme was also well-balanced with respect to the lake at rest as well as moving steady states.However, this scheme suffers from two limitations.Only the semi-discrete version of it was shown to be energy stable.Furthermore, the multi-dimensional version of the scheme was restricted to Cartesian grids.
The main aim of this paper is to design a scheme for the shallow-water equations with bottom topography that is • Fully discrete.
• Able to handle unstructured grids (in several space dimensions).
• Well-balanced with respect to the lake at rest steady state.
To this end, we will design a shock-capturing space-time discontinuous Galerkin (DG) method to approximate the shallow-water equations with bottom topography.The scheme will be an extension of a space-time DG method for conservation laws, proposed in [6,7] and based on earlier work [8] and references therein.The main novelty of the paper will be to handle the bottom topography source terms such that the resulting scheme continues to be energy stable and is in addition, well-balanced.

2.
The space-time DG formulation.The aim of this section is to present the entropy stable space-time DG method for the shallow water equations with variable bottom topography.We start with the description of the mesh.
2.1.The mesh.At the n-th time level t n , we denote the time step as ∆t n = t n+1 − t n and the update time interval as I n = [t n , t n+1 ).For simplicity, we assume that the spatial domain Ω ⊂ R d is polyhedral and divide it into a triangulation T , i.e. a set of open convex polyhedra K ⊂ R d with plane faces.Furthermore, we assume mesh regularity [8] and quasiuniformity.For a generic element (cell) K, we denote ∆x K = diam(K) (element width), The mesh width of the triangulation is ∆x(T ) = max K ∆x K .A generic space-time element is the prism: We also assume that there exists an (arbitrarily large) constant C > 0 such that for all time levels n.

2.2.
The variational formulation.The formulation reads: find V ∆x ∈ V p such that for all W ∆x ∈ V p , where V p denotes the space of piecewise polynomials associated with the space-time mesh with total degree at most p.The variational formulations consists of three quasilinear forms, which we will describe in the following.

2.3.
The DG quasilinear form.The form B DG is given by where the sums are taken over the time steps n = 0, . . ., N −1 and all the cells K ∈ T in the mesh.The last sum is also taken over all neighbouring cells K ′ ∈ N (K) of the cell K.In addition, we have employed the notation The additional boundary terms for the shallow water equations with bottom topography are given by , where for any quantity q, q = 1 2 (q − + q + ) and [[q]] = q + − q − denotes the mean and the jump, respectively.
Upwind fluxes are used for the temporal numerical fluxes This ensures causality and further allows to perform marching in time.
We use an entropy-stable numerical flux given by , which consists of an entropy-conservative flux and a diffusion operator.
Note that the diffusion is added in terms of entropy variables.This is important both for the entropy stability and the well-balancedness of the scheme.
More precisely, we will use a Rusanov type of diffusion, which is given by , where λ max (U; ν) is the maximal wave speed in the direction of ν.
An entropy-conservative flux for the shallow water equations is given by (see [4]) Inserting the numerical fluxes into the form B DG (7), we obtain 2.4.Streamline diffusion operator.Following [7], we need a streamline diffusion operator to minimize oscillations that might arise in the pure space-time DG formulation.However, the streamline diffusion operator of [7] needs to be adapted to balance laws.The equation residual (or intra-element residual) is now The streamline diffusion operator is then given by where L(V ∆x , W ∆x ) denotes a linearised form of the equation in the following sense: For conservation laws this operator can be chosen as The difficulty for balance laws is to include the source term in an appropriate manner without destroying the linearity.For the shallow water equations with bottom topography, it can be chosen as This is a consequence of the following.The residual satisfies and we have This implies and thus, L satisfies the assumptions.The scaling matrix in the streamline diffusion operator remains unchanged.Here, C SD is a positive constant.
2.5.Shock-capturing operator.As in [7], we further need a shock-capturing operator to stabilize any possible oscillations.The shock-capturing operator of [7] need to be modified for the shallow water equations.It is given by with U V = U V ( V n,K ) for brevity and being the cell average.The scaling factor is chosen as 1), and C SC and CSC being positive constants.The scaling factor relies on the integrated intra-element residual Res n,K := and on the integrated boundary residual which needed to be adapted at the spatial boundaries.
3. Entropy stability.We will show that the scheme is entropy-stable even in the presence of source terms due to the bottom topography.
Theorem 3.1.Consider the space-time DG scheme (6) for the shallow water equations with non-constant bottom topography (1).For simplicity, assume that the exact and approximate solutions have compact support inside the spatial domain Ω.Then, the resulting scheme is entropy-stable, i.e. the approximate solutions Proof.Let us use to denote the quantities used without bottom topography.In this case the entropy is given by S(U) = 1 2 hu 2 + 1 2 gh 2 , with the entropy variables being From the entropy stability proof we know (compare [7,6]) We will show This implies The difference in the spatial interior part is given by The difference in the spatial boundary part without diffusion is given by For the diffusion part we can proceed as in the case of the conservation law to obtain The difference in the temporal interior part is given by The difference in the temporal boundary part is given by Combining the equations ( 17), ( 18), ( 19), (20), and (21), we obtain (15).For the streamline diffusion and the shock-capturing term, we can proceed anlogously as for conservation laws to obtain and Combining ( 6), ( 16), (22), and (23), we obtain the stability estimate (14).
4. Well-balancedness.We will show that the space-time DG scheme is wellbalanced.
Theorem 4.1.Consider the space-time DG scheme (6) for the shallow water equations (1).The resulting scheme is well-balanced, i.e. it can preserve the "lake at rest" state exactly.
Proof.We will concentrate on the one-dimensional case, as the proof can be extended to the two-dimensional case in a straightforward manner.Let us consider initial conditions, such that the "lake at rest" condition is satisfied: We consider and we will show that B(V ∆x , W ∆x ) = 0 is satisfied for all W ∆x ∈ V p .
We will use the notation U ∆x = U(V ∆x ) = (h, hu) T for simplicity.Let us start with the temporal part of the DG operator Next, we consider the spatial parts of the DG operator.We need to integrate by parts We start with the resulting interior part and observe that by ( 24).This trivially implies Now we consider the boundary fluxes without the diffusion part and obtain

Let us add and subtract h[[h]] to obtain
Now we consider the diffusion part of the flux: This yields and Combining ( 25), ( 28), (30), and (32) we obtain that the DG form is satisfied, i.e.
For the streamline and shock-capturing terms, first note that the interior residual vanishes by ( 27): Similarly, also the boundary residual vanishes by ( 29) and (31): This trivally implies We conclude that the scheme is indeed well-balanced.
Remark 1.The proof relies on an integration by parts in the spatial terms.If numerical quadrature is used, this introduces an additional error and the resulting scheme will in general not be well-balanced.Therefore, rather than implementing (8) directly, one should implement it after an integration by parts is performed as in (26).This will guarantee that the scheme is well-balanced even if numerical quadrature is used.
5. Numerical Experiments.We present several numerical experiments to demonstrate the shock-capturing streamline diffusion space-time DG scheme.In all experiments, we have chosen to represent the bottom topography exactly.The constant in the scheme are set to C SD = 10, C SC = 1 and CSC = 0.A value of g = 9.812 is used for the gravitational constant in all the experiments.5.1.Lake at rest.Let us consider the following problem taken from [4].The bottom topography is given The initial condition are "lake at rest", i.e. u = 0 and h + b = 1.We compute on the domain [0, 20] and calculate up to the final time t = 10.The bottom topography and the computed surface height of one example run are shown in Figure 1.
The well-balanced property of the scheme is best demonstrated by tabulating relative L 1 errors of the total height h + b (Table 1) with respect to the resolution as well as the polynomial order.
The results show that scheme preserves the exact steady state up to machine precision.The error is slightly bigger for higher polynomial degrees.This is to be expected as the conditioning becomes slightly worse with a higher polynomial degree (or increased resolution).N c p = 0 p = 1 p = 2 20 1.1e-17 3.3e-16 3.3e-15 40 5.6e-18 1.9e-15 1.6e-15 80 2.8e-18 5.0e-15 5.5e-15 160 4.2e-18 1.2e-14 1.3e-14 320 6.9e-18 2.5e-14 1.6e-14 Table 1.Relative L 1 errors of the water surface h + b for the lake at rest case.5.2.Perturbed lake at rest.The main interest is not in steady state itself (as it is known in advance), but rather waves, representing small perturbations of the lake at rest.Therefore, we perturb the lake of rest initial condition of the last experiment by choosing (as in [4]) The bottom topography and the initial velocity are unchanged, i.e. u = 0 and b is given by (34).The exact solution consists of a left-going and a right-going wave.We evolve the flow only up to t = 1.5 such that the waves have not yet left the domain.Figure 2 shows the water surface height h + b for piecewise constant, linear, and quadratic functions.For all polynomial degrees, the waves can be clearly identified without any spurious numerical artifacts or additional waves.The waves are quite diffused in the case of p = 0.However, the accuracy improves to a great extent by considering higher polynomial degrees.We compute up to the final time t = 1.We start on a triangular mesh with 178 cells and perform uniform refinements.Table 2 shows the resulting relative L 1    2. Relative L 1 errors of the water surface h + b for the two-dimensional lake at rest case.reflections.The rightgoing wave is slowed down in the middle by the bottom topography and this then creates a complicated pattern.Figure 4 shows the flow at various times for piecewise constant and piecewise linear functions.The results for piecewise quadratic functions can be found in Figure 5.The results for p = 0 are very smeared, but still qualitatively the flow is captured.The results for higher degrees are much more accurate, even on a coarse mesh.This experiment clearly demonstrates the ability of the proposed space-time DG method to compute small perturbations of the lake at rest steady state, in a stable and accurate manner.

Conclusion.
We extend the shock-capturing streamline diffusion space-time DG method of [6,7] to approximate the shallow water equations with variable bottom topography.The proposed schemes are energy stable, even for the fully discrete case.Furthermore, we show that the proposed schemes are well-balanced i.e. they preserve the lake at rest steady state.The key design elements of the scheme are the use of entropy (energy) variables as degrees of freedom, energy conservative numerical fluxes and discretizations of the bottom topography and the design of suitable streamline diffusion and shock-capturing operators.Numerical experiments that demonstrate the well-balancing of the schemes and the stability and accuracy in resolving small perturbations of the lake at rest steady state are presented.They illustrate that the schemes, particularly with piecewise linear and piecewise quadratic basis functions, perform very well.
The proposed schemes are implicit in time and involve the solution of non-linear algebraic equations in every time step.However, the schemes are unconditionally energy stable.Hence, very large time steps can be chosen.This flexibility is particularly advantageous while dealing with problems in oceanography, where gravity waves can lead to excessively small time steps for explicit methods.The schemes uses unstructured grids in two space dimensions.This is a marked advantage, when dealing with flows in domains with complicated boundaries such as shorelines.Furthermore, space-time DG schemes lend themselves to adaptive algorithms (particularly goal oriented adaptive algorithms) in a natural manner, see [6].Hence, these schemes could form an attractive alternative to standard well-balanced schemes for shallow water flows in complex domains that involve multiple time scales.Such extensions will be considered in the future.An entropy pair is given by S(U) = 1 2 hu 2 + 1 2 hv 2 + 1 2 gh 2 + ghb, Q 1 (U) = 1 2 hu 3 + 1 2 huv 2 + gh 2 u + ghbu, Q 2 (U) = 1 2 hu 2 v + 1 2 hv 3 + gh 2 v + ghbv, which leads to the entropy variables Entropy-conservative fluxes are given by (see [4]) The additional boundary terms are given by The function G in the linearised operator L in the streamline diffusion operator (11) can be chosen as This can be seen as follows.The residual satisfies and

Figure 1 .
Figure 1.Water level h + b and bottom topography b for the lake at rest for p = 2 and N c = 80.

Figure 2 .
Figure 2. h + b for the perturbed lake at rest.

Figure 3 .
Figure 3. Bottom topography b for the lake at rest case.

6 Figure 4 .
Figure 4. h + b for the perturbed lake at rest.The number of cells is 182272 and 45568 for p = 0 and p = 1, respectively.

6 Figure 5 .
Figure 5. h + b for the perturbed lake at rest with 11392 cells and p = 2.
(V ∆x ))U x k (V ∆x , x) − S(V ∆x ) = (−b x u − b y v) L satisfies the assumptions.