## Abstract

High order methods are becoming increasingly popular in shallow water flow modeling motivated by their high computational efficiency (i.e. the ratio between accuracy and computational cost). In particular, Discontinuous Galerkin (DG) schemes are very well suited for the resolution of the shallow water equations (SWE) and related models, being a competitive alternative to the traditional finite volume schemes. In this work, a novel framework for the construction of DG schemes using augmented Riemann solvers is proposed. Such solvers incorporate the source term at cell interfaces in the definition of the Riemann problem, allowing definition of two different inner states in the so-called star region. The benefits of this family of solvers lie in the exact preservation of the Rankine–Hugoniot condition at cell interfaces at the discrete level, ensuring the preservation of equilibrium solutions (i.e. the well-balanced property) without requiring extra corrections of the numerical fluxes. The semi-discrete DG operator becomes nil automatically under equilibrium conditions, provided the use of suitable quadrature rules. The proposed scheme is applied to the Burgers' equation with geometric source term and to the SWE. The numerical results evidence that the proposed scheme achieves the prescribed convergence rates and preserves the equilibrium states of relevance with machine precision.

## INTRODUCTION

Many engineering and environmental problems involve the study of steady and transient free surface water flows where the vertical scale is much smaller than the horizontal dimensions. Such phenomena are usually described by the shallow water equations (SWE) (Cunge *et al.* 1980), a depth-averaged model that considers a hydrostatic pressure distribution in the vertical direction. In realistic scenarios, the flow is often dominated by the so-called source terms (e.g. bed topography, friction, Coriolis forces, etc.), which may appear as a source of mass and/or momentum on the right-hand side of the equations. The methods herein presented focus on the numerical resolution of the SWE with bed slope source term, which is a non-linear, non-homogeneous hyperbolic system of conservation laws.

Suitable discretization techniques for the source term are thus required to provide physically feasible solutions. For a numerical scheme, the preservation of equilibrium solutions is a feature of paramount importance, as many realistic problems of physical interest can be regarded as perturbations over an equilibrium state (Chertock *et al.* 2017). For the SWE system, the most relevant equilibrium condition is known as the lake-at-rest equilibrium (Bermudez & Vázquez-Cendón 1994). This work aims at the design of a novel scheme for the resolution of the SWE and related models in the framework of the Discontinuous Galerkin (DG) method. To this end, the equilibrium condition is generalized for an arbitrary number of spatial dimensions and equations, so that the methods can be presented in a more general form. As a result, the proposed scheme can also be applied to other problems such as the Burgers' equation with source term.

The DG method is a common choice for the numerical resolution of hyperbolic partial differential equations. The origin of this method is usually attributed to Reed and Hill in a paper published in 1973 on the numerical approximation of the neutron transport equation (Reed & Hill 1973). Later on, Cockburn and Shu presented the fundamentals of the scheme in a series of papers published in the 1980s and the 1990s (Cockburn & Shu 1989; Cockburn & Shu 1990). After some years, the DG method became more popular and was successfully applied to the resolution of evolution problems in continuum mechanics, among others. The DG scheme can be seen as a compromise between finite elements (compactness of the stencils, easy -refinement and parallel scalability) and finite volumes (intrinsic conservativeness, the presence of discontinuities between cells and the use of numerical fluxes). Some of the advantages of the DG method, namely its high accuracy, inherent conservativeness, compactness, high parallel efficiency, flexibility for -adaptivity and possibility of using arbitrary geometry and meshes, make it very well suited for the numerical resolution of the SWE (Xing *et al.* 2010; Caviedes-Voullième *et al.* 2020), which is the focus of this work.

The most relevant characteristic that finite volumes and DG methods share may presumably be the use of numerical fluxes at cell interfaces, which allow to exchange information between cells. A suitable calculation of such fluxes is essential for an accurate approximation of the numerical solution ensuring conservation and equilibrium properties. In a DG and/or finite volume framework, numerical fluxes are usually computed by means of the resolution of the Riemann problem (RP) at cell interfaces, which is defined by the original set of equations together with a piecewise constant initial condition given by the DG reconstructed data at cell interfaces. Under equilibrium, the Rankine–Hugoniot (RH) condition must hold across discontinuities in the solution (e.g. steady shock waves or contact waves). In the presence of geometric source terms, the piecewise DG representation of the geometric variable (i.e. polynomial reconstruction with discontinuities at cell interfaces) involves the presence of contact waves at cell interfaces. In this case, the RH condition states that the difference of fluxes across the interface must balance the contribution of the source term, introduced by such a contact wave (Murillo & García-Navarro 2010; Rosatti & Begnudelli 2010; LeVeque 2011).

Most traditional approaches found in the literature for the exact preservation of the RH condition (and ultimately of the equilibrium solution) in the presence of geometric source terms are based on extra modifications of the numerical fluxes (e.g. using the hydrostatic reconstruction, enforcing continuity of the bed elevation at cell interfaces, etc.), which proved to be very robust and useful (Audusse *et al.* 2004; Xing & Shu 2006). However, the solvers therein used do not explicitly account for the full eigenstructure of the system, i.e. they do not account for the contact wave due to the source term in the numerical solution of the RP. To address this issue, augmented Riemann solvers were introduced in an attempt to represent the effect of the source term in the solution of the RP and automatically fulfill RH conditions. When dealing with geometric source terms, it was shown that the consideration of the source term in the definition of the RP is more convenient in order to ensure certain properties of the numerical solution (George 2008; Murillo & García-Navarro 2010). Augmented solvers have been successfully applied to high order schemes (e.g. WENO-ADER and WENO-RK type schemes) (Navas-Montilla & Murillo 2016; Murillo *et al.* 2019; Navas-Montilla *et al.* 2019), allowing well-balanced and energy-balanced numerical solutions that hold the RH condition at cell interfaces with machine precision under equilibrium without extra corrections. Moreover, the augmented solver approach has also been useful to circumvent numerical shockwave anomalies in the finite volume framework (Navas-Montilla & Murillo 2019) and to model realistic scenarios with complex physics (Franzini & Soares-Frazao 2018; Gordillo *et al.* 2019; Martínez-Aranda *et al.* 2019).

As mentioned before, the preservation of certain equilibrium states at the discrete level is a fundamental requirement for a numerical scheme. When considering the SWE with bottom topography, the relevant equilibrium state is the so-called quiescent equilibrium, also known as ‘lake-at-rest’, where the water surface elevation is constant. Numerical schemes satisfying such an equilibrium state are called well-balanced methods (Bermudez & Vázquez-Cendón 1994; Greenberg & Leroux 1996). There is a large variety of well-balanced methods based on Riemann solvers that ensure the preservation of the still water steady state (García-Navarro & Vázquez-Cendón 2000; Alcrudo & Benkhaldoun 2001; Burguete & García-Navarro 2001, 2004; Sanders *et al.* 2003; Chinnayya *et al.* 2004; Yulong & Shu 2005; LeFloch & Thanh 2007, 2011; Bernetti *et al.* 2008; Catella *et al.* 2008; Kesserwani *et al.* 2008; Lee & Wright 2010; Rosatti & Begnudelli 2010; Haleem *et al.* 2015). In 2015, Caleffi *et al.* present an exhaustive study on the treatment of the bed step source term for well-balanced DG schemes, both in the framework of the classical finite volume approach and path-conservative methods. They examine the technique by Kesserwani & Liang (2011) that modifies the discrete bed profile data imposing continuity at cell interfaces, the widespread hydrostatic reconstruction method (Audusse *et al.* 2004) and a path-conservative scheme based on the Dumbser–Osher–Toro (DOT) Riemann solver (Dumbser & Toro 2011). The latter is based on rewriting the SWE as a non-conservative system which includes the trivial equation for the temporal variation of the bed elevation (i.e. representing the conservative product). This approach has been successfully applied to the SWE and enhanced models (Escalante *et al.* 2019) and represents an alternative to augmented solvers, explicitly including the contact wave of the source term in the eigenstructure of the Jacobian. It must be noted that for a hyperbolic system which includes non-conservative products, the choice of the correct path, necessary to define the correct weak solution, is generally non-trivial when complex equilibrium states are sought (Caleffi *et al.* 2015).

The application of augmented Riemann solvers for the computation of the numerical fluxes at cell interfaces represents an alternative to the aforementioned methods and, to the knowledge of the authors, it has never been used in a DG framework. The work herein described focuses on the construction of a DG scheme using augmented Riemann solvers for the resolution of hyperbolic conservation laws with geometric source terms, motivated by the benefits of such family solvers. It is herein evidenced that DG schemes are automatically well balanced when using augmented solvers (i.e. when the numerical fluxes satisfy the RH condition at cell interfaces), provided (i) the use of an exact quadrature rule for the surface and volume integrals and (ii) the reconstructed data satisfy the discrete equilibrium. The time stepping is carried out using an explicit strong stability preserving (SSP) Runge–Kutta (RK) integrator. The proposed methods are applied to the one-dimensional (1D) Burgers' equation with source term and to the SWE with bed topography. For the latter, the HLLS solver is applied, which approaches the eigenstructure of the system by two non-linear waves associated with the convective flux and an extra contact wave associated with the bed discontinuity at cell interfaces. The performance of the HLLS solver in the context of DG-RK schemes is assessed when considering an additional transport equation in the SWE, thus an additional contact wave.

The analytical derivation of the proposed scheme, as well as the proof of the well-balanced property, is carried out in 1D and extended in the Supplementary material, Appendix B, for a general hyperbolic system of equations without any restriction on the spatial dimension. Furthermore, the concept of isotropic equilibrium is introduced as a generalization of the hydrostatic equilibrium. It appears when the matrix of physical fluxes, evaluated at an equilibrium state, is spherical (i.e. proportional to the identity matrix).

The performance of the proposed schemes is herein assessed using different test cases, including steady and transient solutions, that allow to evaluate the well-balanced property and the empirical rates of convergence. The numerical results evidence that the proposed schemes preserve the quiescent equilibrium with machine precision and achieve the expected convergence rates up to third order of accuracy, limited here by the use of the SSPRK3 method.

The paper is organized as follows. In the ‘Mathematical model’ section, the mathematical formulation of the problems of interest is presented. In the ‘General formulation of the well-balanced DG-RK method’ section, DG schemes are introduced and the proposed DG-SSPRK3 augmented scheme is presented. The ‘Application to the Burgers’ equation’ and ‘Application to the SWE’ sections include the numerical results of the application of the proposed scheme to Burgers' and the SWE, respectively, and finally, in the last section, the concluding remarks are presented.

## MATHEMATICAL MODEL

*z*is the bed elevation in the SWE). In this work, is defined as follows:where is the identity matrix,

*u*is a component of related to the equilibrium as it will be shown later (e.g. the water depth,

*h*, in the SWE) and is a non-linear function of

*u*. We will generally consider that , where

*a*is constant.

## GENERAL FORMULATION OF THE WELL-BALANCED DG-RK METHOD

*N*volume cells, denoted by , such that . The DG method can be obtained by applying a traditional Galerkin projection method to each element. Inside each element, the solution is approximated by a linear combination of basis functions as follows:where is the

*l*-th degree of freedom and is the number of degrees of freedom on each element. It is worth noticing that the basis functions only depend upon space and the degrees of freedom are a function of time only.

*q*; is the Gaussian weights and is the determinant of the Jacobian matrix associated with the change of coordinates from the reference element, , to the real element, . In the 1D cases presented in this work, we consider , provided the reference interval is .

### Augmented Riemann solvers

When dealing with hyperbolic conservation laws with geometric source terms (e.g. Equation (1)), non-homogeneous RPs naturally arise at quadrature points across cell interfaces motivated by the discontinuous representation of the solution provided by the DG approach. In such cases, the numerical fluxes at cell interfaces, and , must be computed using augmented Riemann solvers.

*i*and , respectively; only includes the -component and is a singular source term defined at the interface such that , with the property .

#### Resolution of the augmented RP for a scalar equation

#### Resolution of the augmented RP for a system of conservation laws: the HLLS solver

### Well-balancing the scheme

Well-balanced schemes are able to preserve the isotropic equilibrium in Equation (6) with machine accuracy. This means that the spatial operator in Equation (18) must vanish in such an equilibrium condition. The particular requirements a DG method must fulfill in order to be well-balanced are stated in Xing & Shu (2006) and Caleffi & Valiani (2012) as follows:

Under equilibrium, the reconstructed data (using the DG projection) must satisfy the equilibrium condition: this is achieved by applying the DG reconstruction over the equilibrium variable and the geometric variable, and then computing the conserved variable from the difference between these.

Under equilibrium, the numerical flux at the interface must be equal to the physical flux: this is automatically achieved by augmented solvers.

Under equilibrium, the quadrature rule for the surface and volume integrals must be exact: this is achieved by choosing the adequate number of quadrature points, provided the polynomial order.

Note that other different approaches for the construction of well-balanced DG schemes that do not require some of the conditions above can be found in the literature, in particular Wintermeyer *et al.* (2017) proposed a well-balanced DG method for the SWE with inexact integration. These approaches have recently been extended to solve complex multiphase flow problems (Manzanero *et al.* 2020).

In this section, the aforementioned requirements are proved to be fulfilled. The well-balanced equilibrium condition is first derived using the analytical form of the spatial operator in Equation (18). Then, the conditions for the numerical approximation of the integrals in order to preserve the discrete equilibrium with machine accuracy are stated.

p (degree)
. | Truncation error . | MID (1D) . | MID (2D,3D) . |
---|---|---|---|

1 | 2 | 3 | |

2 | 5 | 6 | |

3 | 8 | 9 |

p (degree)
. | Truncation error . | MID (1D) . | MID (2D,3D) . |
---|---|---|---|

1 | 2 | 3 | |

2 | 5 | 6 | |

3 | 8 | 9 |

### The choice of a suitable orthogonal basis: Legendre polynomials

## APPLICATION TO THE BURGERS’ EQUATION

*e*is constant.

*e*, and the geometric variable,

*z*. Then, the projection of the advected variable (see Equation (9)) is computed as follows:

### Case 1. Preservation of an initial equilibrium

*e*, denoted as , is used. It is defined as an integral error norm as follows:with

*q*looping over the quadrature points.

Results for are presented for the well-balanced and non-well-balanced version of the scheme shown in Table 2. Such a non-well-balanced version of the scheme is constructed by means of a homogeneous flux, that is, without including the source term in the definition of the RP, thus preventing the exact balance between fluxes and sources. It is observed that the well-balanced version of the schemes successfully ensures the well-balanced property with machine accuracy. Differently, the non-well-balanced schemes exhibit substantial errors, which are reduced as the order is increased.

. | Non-well-balanced . | Well-balanced . | ||||
---|---|---|---|---|---|---|

Mesh . | DG2 . | DG3 . | DG4 . | DG2 . | DG3 . | DG4 . |

9.22 × 10^{−03} | 2.78 × 10^{−03} | 7.89 × 10^{−06} | 8.88 × 10^{−16} | 4.44 × 10^{−16} | 4.88 × 10^{−15} | |

2.28 × 10^{−03} | 6.96 × 10^{−04} | 4.90 × 10^{−07} | 1.77 × 10^{−15} | 4.44 × 10^{−16} | 5.77 × 10^{−15} |

. | Non-well-balanced . | Well-balanced . | ||||
---|---|---|---|---|---|---|

Mesh . | DG2 . | DG3 . | DG4 . | DG2 . | DG3 . | DG4 . |

9.22 × 10^{−03} | 2.78 × 10^{−03} | 7.89 × 10^{−06} | 8.88 × 10^{−16} | 4.44 × 10^{−16} | 4.88 × 10^{−15} | |

2.28 × 10^{−03} | 6.96 × 10^{−04} | 4.90 × 10^{−07} | 1.77 × 10^{−15} | 4.44 × 10^{−16} | 5.77 × 10^{−15} |

### Case 2. Perturbation of the equilibrium state

Figure 3 shows the exact solution and numerical results for *e* computed by the different schemes mentioned above. It is observed that the well-balancedness of the scheme is crucial for an adequate resolution of the propagation of the perturbation. The results provided by the non-well-balanced schemes show the presence of spurious waves (which are not noticeable for the fourth order scheme) that are only reduced with mesh refinement and by increasing the order of accuracy. Note that the exact solution displayed in Figure 3 is computed using the procedure in Algorithm 1 (Supplementary material, Appendix A), which is based on a fine evolution of the solution using the conservation property along the characteristic curves.

A representation of the evolution of *e* in time is depicted in Figure 4. The numerical results evidence that two different spurious perturbations appear in the non-well-balanced solutions: (i) a traveling perturbation that moves at velocity *u* and (ii) a steady perturbation over the geometrical hump around . Such spurious waves do not appear in the well-balanced solutions.

### Case 3. Convergence rate test

*et al.*2019). A mesh refinement ratio of 2 is considered, leading to the following set of mesh's number of cells: . For two arbitrary meshes and , with the number of cells, the norm of the difference of the computed

*u*in such meshes is calculated as follows:with ,where stands for the grid in which the solution is computed and

*l*is the cell number. Additionally, the error norm can also be defined as follows:

*u*. The numerical results evidence that the prescribed order of accuracy is achieved.

. | WB DG2-SSPRK3 . | WB DG3-SSPRK3 . | ||||||
---|---|---|---|---|---|---|---|---|

Meshes () . | . | . | . | . | . | . | . | . |

2.12 × 10^{−05} | – | 1.27 × 10^{−04} | – | 2.14 × 10^{−07} | – | 1.08 × 10^{−06} | – | |

2.95 × 10^{−06} | 2.85 | 1.85 × 10^{−05} | 2.78 | 8.80 × 10^{−09} | 4.61 | 5.71 × 10^{−08} | 4.25 | |

4.18 × 10^{−07} | 2.82 | 2.53 × 10^{−06} | 2.86 | 6.01 × 10^{−10} | 3.87 | 4.05 × 10^{−09} | 3.82 | |

6.48 × 10^{−08} | 2.69 | 3.78 × 10^{−07} | 2.74 | 6.06 × 10^{−11} | 3.31 | 4.09 × 10^{−10} | 3.31 | |

1.15 × 10^{−08} | 2.49 | 7.03 × 10^{−08} | 2.42 | 7.16 × 10^{−12} | 3.08 | 4.80 × 10^{−11} | 3.09 | |

2.375 × 10^{−09} | 2.28 | 1.56 × 10^{−08} | 2.16 | 8.84 × 10^{−13} | 3.02 | 5.80 × 10^{−12} | 3.05 | |

5.39375 × 10^{−10} | 2.13 | 3.73 × 10^{−09} | 2.06 | 1.10 × 10^{−13} | 3.01 | 7.12 × 10^{−13} | 3.03 |

. | WB DG2-SSPRK3 . | WB DG3-SSPRK3 . | ||||||
---|---|---|---|---|---|---|---|---|

Meshes () . | . | . | . | . | . | . | . | . |

2.12 × 10^{−05} | – | 1.27 × 10^{−04} | – | 2.14 × 10^{−07} | – | 1.08 × 10^{−06} | – | |

2.95 × 10^{−06} | 2.85 | 1.85 × 10^{−05} | 2.78 | 8.80 × 10^{−09} | 4.61 | 5.71 × 10^{−08} | 4.25 | |

4.18 × 10^{−07} | 2.82 | 2.53 × 10^{−06} | 2.86 | 6.01 × 10^{−10} | 3.87 | 4.05 × 10^{−09} | 3.82 | |

6.48 × 10^{−08} | 2.69 | 3.78 × 10^{−07} | 2.74 | 6.06 × 10^{−11} | 3.31 | 4.09 × 10^{−10} | 3.31 | |

1.15 × 10^{−08} | 2.49 | 7.03 × 10^{−08} | 2.42 | 7.16 × 10^{−12} | 3.08 | 4.80 × 10^{−11} | 3.09 | |

2.375 × 10^{−09} | 2.28 | 1.56 × 10^{−08} | 2.16 | 8.84 × 10^{−13} | 3.02 | 5.80 × 10^{−12} | 3.05 | |

5.39375 × 10^{−10} | 2.13 | 3.73 × 10^{−09} | 2.06 | 1.10 × 10^{−13} | 3.01 | 7.12 × 10^{−13} | 3.03 |

## APPLICATION TO THE SWE

*h*, the unit discharge in

*x*, , and the unit discharge in the transverse direction, ; is the physical flux;

*u*is the depth-averaged velocity and m/s

^{2}is the acceleration of gravity. The unit discharge, , is also denoted by

*q*. The source term accounts for the variations in bed geometry:where

*z*is the geometric variable, representing the bed elevation. Note that the third equation in the system in (52) involves a contact wave in the eigenstructure of the Jacobian matrix of and is associated with the passive transport of with celerity

*u*.

Recall that the conservation of the specific discharge in the momentum equation is affected by the use of different left and right numerical fluxes at cell interfaces, allowing to ensure the exact equilibrium between momentum fluxes and the force exerted by the change of bathymetry. On the other hand, the conservation of mass, i.e. variable *h*, is not affected as the source term only acts in the momentum equation.

### Case 1. Preservation of an initial equilibrium

*z*is defined as follows:is considered in the domain and is depicted in Figure 5. The solution is computed at using the DG2 and DG3 schemes in combination with an SSPRK3 integrator, setting and and . For the latter mesh, more than 2,000 time steps are taken. The initial condition is in equilibrium and must thus be preserved by the well-balanced version of the schemes. To assess the performance of the schemes, the infinity error norm for the equilibrium variable, , is used.

Results for are presented for the well-balanced version of the schemes shown in Table 4. It is observed that the proposed schemes successfully ensure the well-balanced property with machine accuracy.

Mesh . | DG2-SSPRK3 . | DG3-SSPRK3 . |
---|---|---|

2.22 × 10^{−16} | 1.00 × 10^{−16} | |

8.88 × 10^{−16} | 3.11 × 10^{−15} |

Mesh . | DG2-SSPRK3 . | DG3-SSPRK3 . |
---|---|---|

2.22 × 10^{−16} | 1.00 × 10^{−16} | |

8.88 × 10^{−16} | 3.11 × 10^{−15} |

### Case 2. Convergence rate test

*h*and , involved in the non-linear dynamics of Equation (52). To empirically assess the order of convergence of the proposed scheme, the following initial condition is chosen:which is defined in the domain . The topography is given by the smooth function . The solution is computed at () using the HLLS, DG2 and DG3 schemes in combination with an SSPRK3 integrator, in eight different meshes with the number of cells 20, 40, 80, 160, 320, 640, 1,280, 2,560. The numerical solution using 20 cells is depicted in Figure 6. Tables 5 and 6 show the relative numerical errors, computed using the and error norms, and empirical convergence rates for

*h*and , respectively. The numerical results evidence that the prescribed order of accuracy is achieved.

. | . | . | ||||||
---|---|---|---|---|---|---|---|---|

Meshes () . | . | . | . | . | . | . | . | . |

1.47 × 10^{−05} | – | 7.17 × 10^{−05} | – | 6.30 × 10^{−05} | – | 4.22 × 10^{−04} | – | |

7.98 × 10^{−07} | 4.20 | 5.34 × 10^{−06} | 3.74 | 3.54 × 10^{−06} | 4.15 | 3.18 × 10^{−05} | 3.73 | |

4.67 × 10^{−08} | 4.09 | 3.47 × 10^{−07} | 3.94 | 2.09 × 10^{−07} | 4.08 | 1.66 × 10^{−06} | 4.25 | |

4.17 × 10^{−09} | 3.48 | 2.97 × 10^{−08} | 3.54 | 1.86 × 10^{−08} | 3.48 | 1.40 × 10^{−07} | 3.571 | |

4.72 × 10^{−10} | 3.14 | 3.30 × 10^{−09} | 3.17 | 2.11 × 10^{−09} | 3.14 | 1.55 × 10^{−08} | 3.17 | |

5.74 × 10^{−11} | 3.03 | 4.00 × 10^{−10} | 3.04 | 2.57 × 10^{−10} | 3.03 | 1.88 × 10^{−09} | 3.04 | |

7.12 × 10^{−12} | 3.01 | 4.93 × 10^{−11} | 3.02 | 3.19 × 10^{−11} | 3.01 | 2.33 × 10^{−10} | 3.01 |

. | . | . | ||||||
---|---|---|---|---|---|---|---|---|

Meshes () . | . | . | . | . | . | . | . | . |

1.47 × 10^{−05} | – | 7.17 × 10^{−05} | – | 6.30 × 10^{−05} | – | 4.22 × 10^{−04} | – | |

7.98 × 10^{−07} | 4.20 | 5.34 × 10^{−06} | 3.74 | 3.54 × 10^{−06} | 4.15 | 3.18 × 10^{−05} | 3.73 | |

4.67 × 10^{−08} | 4.09 | 3.47 × 10^{−07} | 3.94 | 2.09 × 10^{−07} | 4.08 | 1.66 × 10^{−06} | 4.25 | |

4.17 × 10^{−09} | 3.48 | 2.97 × 10^{−08} | 3.54 | 1.86 × 10^{−08} | 3.48 | 1.40 × 10^{−07} | 3.571 | |

4.72 × 10^{−10} | 3.14 | 3.30 × 10^{−09} | 3.17 | 2.11 × 10^{−09} | 3.14 | 1.55 × 10^{−08} | 3.17 | |

5.74 × 10^{−11} | 3.03 | 4.00 × 10^{−10} | 3.04 | 2.57 × 10^{−10} | 3.03 | 1.88 × 10^{−09} | 3.04 | |

7.12 × 10^{−12} | 3.01 | 4.93 × 10^{−11} | 3.02 | 3.19 × 10^{−11} | 3.01 | 2.33 × 10^{−10} | 3.01 |

. | . | . | ||||||
---|---|---|---|---|---|---|---|---|

Meshes () . | . | . | . | . | . | . | . | . |

4.43 × 10^{−04} | – | 2.30 × 10^{−03} | – | 2.20 × 10^{−03} | – | 1.16 × 10^{−02} | – | |

7.80 × 10^{−05} | 2.51 | 4.01 × 10^{−04} | 2.49 | 3.63 × 10^{−04} | 2.57 | 2.20 × 10^{−03} | 2.39 | |

1.08 × 10^{−05} | 2.84 | 5.48 × 10^{−05} | 2.87 | 5.04 × 10^{−05} | 2.85 | 3.31 × 10^{−04} | 2.74 | |

1.57 × 10^{−06} | 2.79 | 7.31 × 10^{−06} | 2.90 | 7.31 × 10^{−06} | 2.79 | 5.16 × 10^{−05} | 2.68 | |

2.68 × 10^{−07} | 2.56 | 1.13 × 10^{−06} | 2.69 | 1.22 × 10^{−06} | 2.57 | 8.73 × 10^{−06} | 2.56 | |

5.51 × 10^{−08} | 2.23 | 2.24 × 10^{−07} | 2.34 | 2.47 × 10^{−07} | 2.31 | 1.74 × 10^{−06} | 2.32 | |

1.27 × 10^{−08} | 2.12 | 5.09 × 10^{−08} | 2.14 | 5.65 × 10^{−08} | 2.13 | 3.98 × 10^{−07} | 2.13 |

. | . | . | ||||||
---|---|---|---|---|---|---|---|---|

Meshes () . | . | . | . | . | . | . | . | . |

4.43 × 10^{−04} | – | 2.30 × 10^{−03} | – | 2.20 × 10^{−03} | – | 1.16 × 10^{−02} | – | |

7.80 × 10^{−05} | 2.51 | 4.01 × 10^{−04} | 2.49 | 3.63 × 10^{−04} | 2.57 | 2.20 × 10^{−03} | 2.39 | |

1.08 × 10^{−05} | 2.84 | 5.48 × 10^{−05} | 2.87 | 5.04 × 10^{−05} | 2.85 | 3.31 × 10^{−04} | 2.74 | |

1.57 × 10^{−06} | 2.79 | 7.31 × 10^{−06} | 2.90 | 7.31 × 10^{−06} | 2.79 | 5.16 × 10^{−05} | 2.68 | |

2.68 × 10^{−07} | 2.56 | 1.13 × 10^{−06} | 2.69 | 1.22 × 10^{−06} | 2.57 | 8.73 × 10^{−06} | 2.56 | |

5.51 × 10^{−08} | 2.23 | 2.24 × 10^{−07} | 2.34 | 2.47 × 10^{−07} | 2.31 | 1.74 × 10^{−06} | 2.32 | |

1.27 × 10^{−08} | 2.12 | 5.09 × 10^{−08} | 2.14 | 5.65 × 10^{−08} | 2.13 | 3.98 × 10^{−07} | 2.13 |

### Case 3. Riemann problems

The numerical resolution of different RPs, provided in Table 7, is considered in this section. The spatial domain is given by , the discontinuity is positioned at and the acceleration of gravity is set to = 9.81 m/s^{2}. Numerical solutions are plotted for , and and compared with exact solutions. For all the problems, we set .

Label . | . | . | . | . | . | . | . | . | T
. |
---|---|---|---|---|---|---|---|---|---|

RP1 | 1.0 | 0.1614067989 | 0.0 | 0.0 | 1.0 | 0.1614067989 | 0.0 | 0.05 | 0.02 |

RP2 | 1.0 | 0.6245627691 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.3 | 0.01 |

RP3 | 1.0 | 1.0 | 0.0 | 0.0 | 1.0 | 2.0 | 0.0 | 0.0 | |

RP4 | 0.4 | 1.6 | 2.3 | 1.686811595600704 | 1.0 | 0.0 | 0.0 | 0.0 | 0.1835 |

Label . | . | . | . | . | . | . | . | . | T
. |
---|---|---|---|---|---|---|---|---|---|

RP1 | 1.0 | 0.1614067989 | 0.0 | 0.0 | 1.0 | 0.1614067989 | 0.0 | 0.05 | 0.02 |

RP2 | 1.0 | 0.6245627691 | 1.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.3 | 0.01 |

RP3 | 1.0 | 1.0 | 0.0 | 0.0 | 1.0 | 2.0 | 0.0 | 0.0 | |

RP4 | 0.4 | 1.6 | 2.3 | 1.686811595600704 | 1.0 | 0.0 | 0.0 | 0.0 | 0.1835 |

Units in m, m^{2}/s and s.

The first RP (RP1) is computed using two different grids, composed of 160 and 320 cells, respectively, to show the convergence to the exact solution with mesh refinement. The solution for the water surface elevation, normal discharge and shear discharge are plotted in Figure 7. It is observed that the intermediate states are accurately predicted, as well as the location of the shock wave and the rarefaction wave. As no slope limiter is used, the second order scheme exhibits a spurious spike in the cell containing the right-moving shock. However, the presence of this spike does not affect the rest of the solution. When moving from the first order scheme to the second order scheme, a significant gain in accuracy is evidenced. A detail of the solution around the contact discontinuity is also presented in Figure 7, showing that the second order scheme accurately reproduces such a discontinuous wave without oscillations, even when the slope limiter is not used.

The second RP (RP2) consists of an equilibrium RP, as the initial data is a discontinuous solution which satisfies the RH condition. The numerical solution, depicted in Figure 8, is computed using 160 cells. As expected, both the first order and second order schemes preserve the steady moving water equilibrium with machine precision.

The fourth RP (RP4) involves the computation of a slowly moving hydraulic jump (strong shock), usually linked to the slowly moving shock anomaly in the framework of FV schemes (Navas-Montilla & Murillo 2019). The aim of this case is to assess the robustness of the augmented DG-RK scheme without the use of a limiter and to determine whether or not this scheme suffers from the slowly moving shock anomaly. Figure 10 shows the numerical solution for *h* and computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 80 cells. It is observed that both schemes produce a spurious oscillation downstream the hydraulic jump, however, of a different nature. In Figure 11, the evolution in time of the average value of provided by the schemes at a cell on the left side of the initial discontinuity () is depicted. It is observed that the first order scheme produces a spurious spike of the unit discharge (up to m^{2}/s) when the shock crosses the aforementioned cell, which corresponds to the slowly moving shock anomaly (Navas-Montilla & Murillo 2019). Contrarily, the DG2-SSPRK3 does not show this anomaly, being the oscillations only produced by the lack of a limiter. This evidences that the proposed DG2-SSPRK3 scheme is free from the slowly moving shock anomaly.

The reason for this is that FV schemes inherently enforce the width of shock waves due to the cell-averaged representation of data, whereas DG schemes do not. When cell-averaging the data around the shock, non-physical values of the conserved variables may appear in the case of a non-monotone (non-linear) Hugoniot locus (Navas-Montilla & Murillo 2019), which eventually produce the observed spike in the SWE. Contrarily, DG methods are based on a modal/nodal reconstruction of pointwise data, without averaging, and therefore do not exhibit this problem. Nevertheless, this behavior should be further investigated in future works.

## CONCLUDING REMARKS

In this work, the application of augmented Riemann solvers in the framework of the DG method is addressed. An explicit DG-SSPRK3 scheme based on augmented Riemann solvers is constructed. The proposed scheme is applied to the resolution of hyperbolic conservation laws with geometric source terms, in particular, to the Burgers' equation with geometric source term and to the SWE with bed topography.

The use of augmented Riemann solvers (e.g. the HLLS solver) proves to be very useful and effective for the preservation of equilibrium states. Thanks to the definition of the RP considering the contribution of the source term at cell interfaces, the exact RH condition is satisfied at cell interfaces by the Riemann solver without requiring extra flux corrections. It is proved that the DG spatial operator becomes nil under equilibrium, provided (i) the use of an exact quadrature rule for the surface and volume integrals and (ii) the reconstructed data satisfy the discrete equilibrium. Such a discrete equilibrium property of the scheme is proved with the independence of the number of spatial dimensions.

For a proper analysis of the fundamental solutions obtained by the combination of augmented solvers with the DG method, the scheme herein presented includes neither extra flux corrections nor slope limiting algorithms. Nevertheless, the incorporation of additional numerical fixes (i.e. slope limiters, wet-dry treatment, etc.) to address realistic problems of engineering and environmental interest is straightforward. Furthermore, the application of this scheme to higher spatial dimensions is also straightforward thanks to the rotational invariance property of the systems (i.e. 1D RPs are always solved, with the independence of the number of spatial dimensions).

The numerical results evidence that the proposed scheme achieves the prescribed convergence rates and preserves the equilibrium states of relevance with machine precision, for the problems considered in this work. Convergence rates have been assessed up to third order of accuracy, as this is the limiting order due to the SSPRK3 time integrator. Arbitrary order numerical schemes can be constructed using the DG-ADER approach (Dumbser & Munz 2005; Fambri *et al.* 2017), where augmented solvers can be applied as presented here for the resolution of the Derivative RP using the LFS solver in Navas-Montilla & Murillo (2016). The results also evidence that the proposed DG2-SSPRK3 scheme is free from the slowly moving shock anomaly.

## ACKNOWLEDGEMENTS

This work was funded by the Spanish Ministry of Science and Innovation under the research project PGC2018-094341-B-I00. This work has also been partially funded by Gobierno de Aragón through Fondo Social Europeo (Feder 2014-2020 ‘Construyendo Europa desde Aragón’).

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2020.206.