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
is the vector of conserved quantities that takes values on 
, with 
is the number of equations. 
is the matrix of physical fluxes. The source term is defined as a geometric source term:
a geometric function that depends upon the position (e.g. z is the bed elevation in the SWE). In this work, 
is defined as follows:
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.
be a steady solution of (1). The equilibrium state will be called isotropic if the physical flux matrix at equilibrium is of the form:
the equilibrium value of one of the components of vector 
that participates in the equilibrium condition. Note that the hydrostatic equilibrium is a particular case of the isotropic equilibrium.
) which participates in the equilibrium condition.
is reduced to vector 
, which stands for the physical flux in the x-direction according to the standard notation and 
will include a single derivative of the geometric variable, 
.GENERAL FORMULATION OF THE WELL-BALANCED DG-RK METHOD
, is discretized in 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:
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.
, and integrating by parts the second term of (10) using Green's first identity:
and 
is bi-valued, the concept of numerical flux has to be used and the previous identity is rewritten as follows:
and 
are the cell entering and leaving numerical fluxes, respectively. Note that we will consider different left and right numerical fluxes computed from the solution of the RP at the interface, denoted by 
and 
, respectively. Using the equality in Equation (12), we can rewrite Equation (10) as follows:
.
, the following property is satisfied
, where 
is the Kronecker delta and 
is the value of the integral over 
. Using the result in (15), (14) becomes a decoupled system of equations for the unknown degrees of freedom 
that reads
. Equation (16) can be rewritten as follows:
, the volume integrals in (18) are computed as follows:
is the number of the Gaussian quadrature points; 
, 
and 
are the evaluation of the flux, source term and the basis function at the Gaussian quadrature points 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 
.
using a TVD RK method. The third order SSP RK method will be herein used:
degrees of freedom at 
, denoted by 
.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.
and 
are the left and right reconstruction of the variables at cell interfaces, provided by the DG polynomials in cells i and 
, respectively; 
only includes the 
-component and is a singular source term defined at the interface such that 
, with the property 
.
), which originates a jump of the fluxes and conserved variables across the interface. The approximate solution, hereafter referred to as 
, is therefore composed by at least (depending on the number of waves) two different internal states in the so-called star region, separated by a contact discontinuity at 
produced by the source term. In the vicinity of 
, left and right states of the approximate solution will be hereafter denoted by 
and 
, respectively, expressed as follows:
can be also constructed from flux values 
and 
with a similar structure than 
. In this case, also intercell values for the fluxes can be defined at both sides of the 
-axis as follows:
and 
), Equation (23) is recovered. Note that the numerical fluxes will be generally expressed using the subscript 
notation as 
and 
.Resolution of the augmented RP for a scalar equation
and
is the approximate wave celerity at the interface.Resolution of the augmented RP for a system of conservation laws: the HLLS solver
waves (e.g. 1D SWE), the HLLS solver accounts for the full eigenstructure of the system. The numerical fluxes are computed as follows:
and 
are the left- and right-hand limits to the cell edge of the conserved quantities (
) and their 
-th time derivatives and 
is a matrix related to the source term. Further details can be found in Navas-Montilla & Murillo (2016).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.
-th cell-leaving numerical fluxes in Equation (18) are equal to the physical fluxes in the cell 
, according to Equation (27). This allows Equation (18) to be rewritten as follows:
is a steady-state equilibrium solution.
Polynomial degree, truncation order and maximum integrand degree (MID) in the spatial operator for 1D and 2D/3D problems
| 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
norm, which satisfies Equation (15). A recursive formula for Legendre polynomials, hereafter denoted by 
, is given by the following equation:
, can be obtained as follows:
APPLICATION TO THE BURGERS’ EQUATION
is the problem variable and 
is a geometric variable. The equilibrium condition for Equation (39) is obtained when 
vanishes, yielding
. Steady solutions of (39) will satisfy that e is constant.
and the approximate wave celerity yields 
. The approximation of the integral of the source term at cell interfaces is selected as follows:
. This allows the numerical fluxes to be written as follows:
Case 1. Preservation of an initial equilibrium
and is depicted in Figure 1. The solution is computed at 
using the DG2, DG3 and DG4 schemes in combination with an SSPRK3 integrator, setting 
and 
and 
. 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 in the preservation of the initial equilibrium, the infinity error norm for the equilibrium variable e, denoted as 
, is used. It is defined as an integral error norm as follows:
quadrature points.
Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
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.
error at 
for the well-balanced and non-well-balanced DG2-SSPRK3, DG3-SSPRK3 and DG4-SSPRK3 schemes in two grids with 
and 
| 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
Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
is given by
. The computational domain is given by 
. The solution is computed at 
(
) using the DG2, DG3 and DG4 schemes in combination with an SSPRK3 integrator, in two different meshes with 
and 
.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.
Numerical solution for e computed by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (purple), DG3-SSPRK3 (yellow) and DG4-SSPRK3 schemes (green) at 
in a 
(left) and 
grid (right). The initial condition is depicted using the black dashed line and the exact solution using the black solid line. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
Numerical solution for e computed by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (purple), DG3-SSPRK3 (yellow) and DG4-SSPRK3 schemes (green) at in a (left) and grid (right). The initial condition is depicted using the black dashed line and the exact solution using the black solid line. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
Evolution of the computed e provided by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (left), DG3-SSPRK3 schemes (right) at 
in a 
grid.
Evolution of the computed e provided by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (left), DG3-SSPRK3 schemes (right) at in a grid.
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
. 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:
,
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:
. The solution is computed at 
(
) using the DG2, DG3 schemes in combination with an SSPRK3 integrator, in the meshes indicated above. Table 3 shows the relative numerical errors, computed using the 
and 
error norms, and empirical convergence rates for u. The numerical results evidence that the prescribed order of accuracy is achieved.Convergence rates for the 
and 
error norms, computed by the well-balanced DG2-SSPRK3 and DG3-SSPRK3 schemes
| 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
-split SWE represent the projection of the 2D SWE system onto the 
-axis, as a result of the rotational invariance property of the SWE (Murillo & Navas-Montilla 2016). They are usually solved at cell interfaces when addressing the computation of 2D shallow water flows. The 1D, 
-split SWE over the irregular bed are expressed as follows:
is the vector of conserved variables, namely the water depth, 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/s2 is the acceleration of gravity. The unit discharge, 
, is also denoted by q. The source term 
accounts for the variations in bed geometry:
and is associated with the passive transport of 
with celerity u.
, the geometric variable, z, and discharge, 
. The degrees of freedom of the water depth will be computed as follows:
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
, is set, where z is defined as follows:
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.
Initial 
(green) and bed elevation z (black).
Numerical solution for 
(left) and 
(right) computed by the well-balanced first order Godunov's scheme (green), DG2-SSPRK3 (purple) and DG3-SSPRK3 (orange) at 
in a 
grid. The reference solution, computed by a third order WENO-ADER scheme (Navas-Montilla & Murillo 2016) in a very fine grid (
), is depicted in black. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
Numerical solution for (left) and (right) computed by the well-balanced first order Godunov's scheme (green), DG2-SSPRK3 (purple) and DG3-SSPRK3 (orange) at in a grid. The reference solution, computed by a third order WENO-ADER scheme (Navas-Montilla & Murillo 2016) in a very fine grid (), is depicted in black. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
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.
error at 
for the well-balanced DG2-SSPRK3 and DG3-SSPRK3 schemes in two grids with 
and 
| 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
, 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:
. 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.Convergence rates for the water depth and discharge, h and 
, associated with the 
and 
error norms, computed by the well-balanced DG3-SSPRK3 schemes
 |  | |||||||
|---|---|---|---|---|---|---|---|---|
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 |
Convergence rates for the water depth and discharge, h and 
, associated with the 
and 
error norms, computed by the well-balanced DG2-SSPRK3 schemes
 |  | |||||||
|---|---|---|---|---|---|---|---|---|
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/s2. Numerical solutions are plotted for 
, 
and 
and compared with exact solutions. For all the problems, we set 
.
Initial condition for the RPs
| 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, m2/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.
RP1. Numerical solution for 
(top), 
(middle) and 
(bottom) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at 
s using 160 (left) and 320 cells (right), compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
RP1. Numerical solution for (top), (middle) and (bottom) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 160 (left) and 320 cells (right), compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
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.
RP2. Numerical solution for 
(left) and 
(right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at 
s using 160 cells, compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
RP2. Numerical solution for (left) and (right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 160 cells, compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
-direction, but there is velocity in the transverse (shear) direction with a discontinuity at the interface. Ideally, the left and right regions of the flow should be decoupled as there is no viscosity in the equations and thus no explicit shear stresses. If a complete Riemann solver such as the (A)Roe or HLLC(S) solver were used, the exact solution would be obtained and the left and right regions would be decoupled in the numerical solution. When using the HLLS solver, which does not explicitly account for the contact wave, the shear discontinuity in 
is smeared out due to the numerical diffusion intrinsic to this particular solver (Navas-Montilla & Juez 2020). The numerical solution for this problem is computed using 80 and 160 cells at times 
and 
s, and is depicted in Figure 9. Only 
is shown in the figure, as the other variables are constant and have no interest. The exact solution and the analytical solution of the equivalent equation of the first order HLLS solver is also depicted. The latter is given by
is the wave celerity (Navas-Montilla & Juez 2020). Figure 9 shows that the first order scheme adds a large amount of numerical diffusion, which would involve a strong mixing between the left and right regions in the context of the resolution of a shear layer. This diffusion increases as the cell size or the simulation time is increased. On the other hand, the DG2-RK3 scheme reproduces a sharper transition across the contact discontinuity (roughly two-cell wide) and shows no dependence on the simulation time (this is noticed as the DG polynomials show no discontinuities across cell interfaces). This is a clear advantage when compared to other high-order numerical methods, such as WENO-ADER schemes, which show a more smeared profile (Navas-Montilla & Juez 2020). A small overshoot is shown in the DG2-RK3 solution as no slope limiter is used. However, the performance of such a scheme is by far better than the first order version of the HLLS solver and even better than other higher order schemes.
RP3. Numerical solution for 
computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at 
(left) and 
s (right) using 80 (top) and 160 cells (bottom), compared with the exact solution (black solid line) and the analytical solution of the equivalent equation for the first order scheme (dashed line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
RP3. Numerical solution for computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at (left) and s (right) using 80 (top) and 160 cells (bottom), compared with the exact solution (black solid line) and the analytical solution of the equivalent equation for the first order scheme (dashed line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
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 
m2/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.
RP4. Numerical solution for h (left) and 
(right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at 
s using 80 cells, compared with the exact solution (black solid line). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
RP4. Numerical solution for h (left) and (right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 80 cells, compared with the exact solution (black solid line). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
RP4. Evolution in time of the average value of 
at a cell on the left side of the initial discontinuity (
) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 scheme (purple). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
RP4. Evolution in time of the average value of at a cell on the left side of the initial discontinuity () computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 scheme (purple). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/hydro.2020.206.
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.
%20cropped.png?versionId=7142)