Well-Balanced High-Order Discontinuous Galerkin Methods for Systems of Balance Laws

This work introduces a general strategy to develop well-balanced high-order Discontinuous Galerkin (DG) numerical schemes for systems of balance laws. The essence of our approach is a local projection step that guarantees the exactly well-balanced character of the resulting numerical method for smooth stationary solutions. The strategy can be adapted to some well-known different time marching DG discretisations. Particularly, in this article, Runge–Kutta DG and ADER DG methods are studied. Additionally, a limiting procedure based on a modified WENO approach is described to deal with the spurious oscillations generated in the presence of non-smooth solutions, keeping the well-balanced properties of the scheme intact. The resulting numerical method is then exactly well-balanced and high-order in space and time for smooth solutions. Finally, some numerical results are depicted using different systems of balance laws to show the performance of the introduced numerical strategy.


Introduction
We consider one-dimensional systems of balance laws of the form ∂ t U + ∂ x F(U) = S(U)H + G(U). (1) Here, U(x, t) is the state vector taking values in the space of admissibly states Ω ⊂ R M ; F : Ω → R M is the flux; S(U) : Ω → R M , is a function that takes into account the geometrical source terms that are multiplied by the known function H : R → R; and G(U) contains algebraic source terms. We suppose that the system is hyperbolic, and therefore the jacobian matrix of the flux has M distinct real eigenvalues.
The system (1) often admits steady-state solutions that, in most cases, are non-trivial. Moreover, these solutions usually convey significant physical meaning. A paradigmatic example of this is the lake-at-rest solution of the well-known shallow-water equations. The correct preservation of these solutions plays an essential role in any simulation. Wellbalanced (or exactly well-balanced) numerical schemes aim to preserve all or at least a relevant set of specific steady-state solutions up to machine precision. A precise definition of well-balanced method is given in Section 2.
An essential advantage of well-balanced schemes is that they can accurately resolve small perturbations to steady-state solutions with relatively coarse meshes. Indeed, if a numerical method cannot preserve stationary solutions or recover a stationary solution after a perturbation, the oscillations introduced may destroy the numerical solution for long integration times. Additionally, the numerical solution may be lost entirely if the perturbation is in the same order as the discretisation errors introduced by a non well-balanced method. This paper is organised as follows. In Section 2, a brief description of the considered numerical schemes are performed, including Runge-Kutta and ADER discretisations in time. Section 3 describes the followed well-balanced methodology proposed for the numerical schemes studied. Next, in Section 4, a limiter procedure that preserves the well-balanced features of the numerical schemes is proposed. Then, in Section 5 some numerical tests are provided for the Shallow water, Burger, Euler with gravity and the hyperbolic dispersive Gavrilyuk-Favrie equations. Finally, some conclusions are given in Section 6.

Discontinuous Galerking Method
We briefly recall the DG method to discretise the general system (1) for one-spatial domains. We consider a set of element cells Ω i = [x i− 1 2 , x i+ 1 2 ], with constant length ∆x. The numerical solution of the PDE (1) at a given time t at each subset Ω i is represented by u h (x, t) and is described in terms of piece-wise polynomials of degree N on the spatial direction. We shall denote by U h the space of piece-wise polynomials up to degree N so that u h (x, t)| Ω i ∈ U h , ∀i.
An orthogonal nodal basis is defined by the Lagrange interpolation polynomials over the (N + 1) Gauss-Legendre quadrature nodes on the cell element Ω i . We stress that the piece-wise solution u h may be discontinuous across elements, allowing jump discontinuities at cell interfaces.
Given the element Ω i , the discrete solution u h projected on the nodal spatial basis polynomials Φ l (x) reads withû i,l (t) being the unknown degrees of freedom, and where the Einstein summation convention over two repeated indices has been considered. In this work, we suppose that H(x) is a given function and it is known at every point of the domain. In the subsequent notation, we omit the space dependency of Φ l (x) for the sake of simplicity. The DG methods results then from multiplying the PDE system (1) with a polynomial test function Φ k ∈ U h and integrating over the spatial element Ω i . That leads to the following weak formulation where we seek u h (x, t) ∈ U h such that, for all Ω i , holds for all test functions Φ k (x) ∈ U h . The dependence on x and t has been dropped for the sake of compactness. Using (2), the flux function is integrated by parts and the weak formulation (3) can be rewritten as where the first integral contains the so-called element mass matrix, which is diagonal since the DG polynomial basis is orthogonal. Likewise, Ω • i denotes the set of all interior points of Ω i . The second integral accounts for the smooth part of the flux, while the last integral accounts for the presence of the continuous source terms. Finally, F i−1/2 − F i+1/2 describes the jump across element interfaces given by a consistent numerical flux of the homogeneous system where u − h,i+1/2 and u + h,i+1/2 are the boundary extrapolated states at the left and right of the interface x i+1/2 , respectively. Note that the numerical fluxes are multiplied by the evaluation of the Lagrange base at the cell interfaces, Φ k (x i±1/2 ). In this work, the standard HLL Riemann solver is used. However, the numerical methods described in this work can be easily extended to more accurate Riemann solvers such as the Polynomial Viscosity Methods or Rational Viscosity Methods (see [86][87][88][89]). Definition 1. Given a continuous stationary solution U of system (1), the semi-discrete numerical scheme (4) is said to be well-balanced (or exactly well-balanced) if the set {û i,l = U(x i,l )} is an equilibrium of the ODE system, where {x i,l } is the set of nodal spatial points. Note that this definition mimics the one given in [35] or [31] for finite volume methods.
It is clear that the semi-discrete method (4), as it is formulated, is not well-balanced in general. The following section will describe a procedure to define DG well-balanced methods.
The semi-discrete numerical scheme (4) can be then approximated with any timemarching numerical method. This work uses classical Runge-Kutta methods and the family of one-step ADER-DG methods based on a predictor-corrector scheme. In the next Subsection, a brief description of the aforementioned time discretisation methods is given. The interested reader can be referred to [38][39][40][41]80,90] for DG Runge-Kutta methods or [67][68][69][70][71][72] for ADER-DG methods.

Runge-Kutta Time Discretization
We first remark that due to the choice of an orthogonal DG polynomial basis, the mass-matrix appearing in the semi-discrete system (4) is diagonal, and that leads us to write (4) for each element Ω i and degree of freedomû i,k in the more compact form where L(û i,k (t)) is the spatial discretization operator. To discretise the time variable t, total variation diminishing (TVD) schemes have been considered in this work. For instance, the third order Runge-Kutta scheme applied to the time interval [t n , t n + ∆t], ∆t being the time-step, readŝ

ADER Time Discretization
The predictor-corrector ADER-DG approach leads to a high-order one-step method in both space and time. Here, the high-order accuracy in time is reached with a local space and time predictor guess. The spatial-time predictor q h (x, t) is computed based on a weak formulation of the PDE. Given the known solution u h (x, t n ) at time t n , a Cauchy problem in the cell element is solved without considering interactions between the neighbour elements. The predictor solution is then used in (4) to compute the numerical solution u h (x, t n+1 ) at the subsequent time step.
Here, we briefly recall the numerical procedure. Given a space-time control volume Ω i × [t n , t n+1 ], the predictor solution q h is projected onto the local space-time basis, where again the Einstein notation has been used. If the system (1) is multiplied by the space-time test function θ k ≡ θ k (x, t) and integrated over Ω i × [t n , t n+1 ], it yields Using Equation (8), Equation (9) becomes a nonlinear system for the unknownsq i l that can be solved via a fixed point algorithm detailed in [65,91]. Notice that (9) is completely local, without any interaction with neighbouring cells. The iterative method converges in N + 1 steps at most for linear homogeneous systems [92].
Finally, the one-step high-order in both space and time ADER-DG method, after replacing u h (x, t) by q h (x, t) in the right-hand side of (4), and integrating in the time interval [t n , t n+1 ], reads where

Well-Balanced Technique
In this Section, we detail a well-balanced technique for the Runge-Kutta DG methods (4)-(7) and the ADER-DG method (9) and (10).

Well-Balanced Runge-Kutta DG Methods
We propose now a spatial DG discretization method that exactly preserves stationary solutions. Later, we will discuss the well-balanced technique for the chosen time marching discretization. As in [35], we begin by computing a suited stationary solution in each element Ω i at each time step t n :

1.
Seek for the stationary solution u * i (x) of the minimization problem in Ω i The previous problem is in general difficult. Nevertheless if u * i (x) depends on a set of parameters, problem (12) reduces, in general, to a non-linear system of equations if the integrals in (12b) are approximated by quadrature formulas defined in terms of the Gaussian nodal points of the DG polynomials. Moreover, in some situations, like the ones considered below, (12a) and (12b) are enough to determine uniquely the local stationary solution. For instance, if u * i (x) is of the form then the integration of (12b) by using quadrature rules leads to where ω l , x l are the Gaussian quadrature rule weights and nodes respectively. This quadrature rule must be of the same order of the DG space approximation at the element Ω i . Similarly, for stationary solutions u * i (x) of the form the constant c i can be simply computed from (12b) and yields 2.
Compute the fluctuation DG polynomial: where To preserve stationary solutions, the numerical method (4) is modified as follows: whereĤ i,l are the coordinates of the projection of H(x) into U h , and is the usual first-order consistent numerical flux evaluated at the states where u * i and u * i+1 are the stationary solutions computed at the step 1 at cells Ω i and Ω i+1 , respectively, and v ) are the extrapolated values of the fluctuation computed at step 2 onto the left and right cell interfaces respectively. Finally, F * i+1/2 denotes the evaluation of the physical flux function at the stationary solution at the intercells:

Remark 2.
The procedure to compute a stationary solution (13)- (16) can be generalized for some other different stationary solutions, or even for a k-parameter family of stationary solutions as described in [35], u * h (x; c i,1 , . . . , c i,k ), c i,r ∈ R.
Remark 3. We stress that, as u * i is a solution of (12), then (18) is of the same order of accuracy than (4) (see in [35] for more details). Theorem 1. The semi-discrete numerical method (18) is well-balanced for a continuous global stationary solution u * , provided that every u * i defined by (12) verifies u * | Ω i = u * i for each Ω i , in the sense thatû * i,k , coordinates of the projection of u * into U h , is an equilibrium of the ODE system (18).
Proof. We consider as initial conditionû i,k (0) =û * i,k , that implies and the three volume integrals are zero. In the same way, the fluctuations v h| Ω i ( , are exactly zero as well. The continuity of u * implies that ) and the consistency of the numerical flux implies that The same reasoning can be applied to prove that F * i+1/2 − F i+1/2 = 0. We conclude thatû * i,k in an equilibrium of the ODE system given by the semi-discrete numerical method (18).
Once the semi-discrete well-balanced method (18) has been proposed, the Cauchy problem for an autonomous ODE system has to be solved. Here L is given by the right-hand side of (18). Clearly, a stationary solution u * h is solution of (21) if and only if Therefore, given a time discretization of the form for (21), the final discrete method is well-balanced for the stationary solution u * h if and only if Finally, it is clear that (22) implies (23) for any Runge-Kutta method as the ones considered in this work.

Well-Balanced ADER DG Methods
To design a well-balanced method based on the ADER time discretization, we follow a similar strategy as for the well-balanced DG Runge-Kutta method. Given an element Ω i , a time-step t n , and the known numerical solution u h (x, t n ) : 1.
First, we seek for the discrete stationary solution u * i (x) of the Cauchy problem in Ω i , as in the first step of the previous Section 3.1.

2.
Then, we define the fluctuation DG polynomial v h and the local space-time predictor given in (9) is modified as follows: Once the nonlinear system for the unknownv i l has been solved via a fixed point algorithm, the space-time predictor q h (x, t) is recovered as

3.
Finally, the one-step ADER-DG method described in (10) is modified as follows: where is the usual first order consistent numerical flux evaluated at the reconstructed predictor states are the stationary solutions computed at the step 1 at cells Ω i and Ω i+1 , respectively, and v − h,i+1/2 (t) and v + h,i+1/2 (t) are the extrapolated values of the fluctuation computed at step 2 onto the left and right cell interfaces, respectively. Additionally, F * i+1/2 denotes the evaluation of the physical flux at the local stationary solution at the inter-cells: . Note that the added terms in the DG system (18) and ADER-DG system (26) are zero up to the order of the method. In this way, when the solution is far from a stationary solution, the solution yielded by these systems will be the same as their non-well-balanced counterparts (4) and (10) up to the order of the method. Additionally, in Section 5 a numerical example is given where a comparison between a well-balanced and non wellbalanced simulation for a problem far from the equilibrium.
Finally, we stress that, although the techniques here have been described for the particular case of one-dimensional problems, all methods can be easily extended in a dimension by dimension fashion for the case of multiple spatial dimensions. For an implementation of the DG method in the higher spatial dimension case, we refer the reader to the work in [76].

A Well-Balanced WENO Limiting Procedure
We now detail the limiting procedure employed in this work, based on Zhong and Shu's work in [83], where the authors propose a weighted essentially non-oscillatory (WENO) limiter for Runge-Kutta DG methods.
The main difference with our proposal is that we limit the fluctuations with respect to the stationary solutions, denoted by v h , rather than the solution itself. Once the fluctuation has been limited, the numerical solution can be recovered by adding the stationary solution.
The limiter can be seen as a two-steps procedure as in [83]: in the first stage, the cells that need limiting, or troubled cells, are identified. Then, in the second step, the cells that were deemed as troubled are limited.
For both Runge-Kutta and ADER time-marching discretisation, the limiting procedure is applied a priori and only once at the beginning of each time step.

Detector of Troubled Cells
We use the TVD minmod limiter detector following ideas detailed in [41,83]. We first definev where v h (x, t n ) denotes the fluctuation with respect to the discrete stationary solution We also define where v + i−1/2 and v − i+1/2 are the evaluation of the polynomial fluctuation v h (x, t n ) onto the left and right cell interfaces respectively. Then, the minmod is applied as follows: where and the following minmod functions are given by Here, M is a parameter to be properly chosen depending on the PDE system. Finally, a cell is marked as troubled if minmod(a, b, c) = a.

Limiting Procedure
The limiting technique for the cells marked as troubled by the detector described in the previous Subsection is now detailed. The general idea is to reconstruct a new polynomial solution at a troubled cell as a convex combination of its neighbouring polynomials.
If the cell Ω i is deemed troubled, then the first step consist on projecting the fluctuation at each neighbour element Ω j , j ∈ {i − 1, i, i + 1}, into a first-order DG polynomial with the following form: that is a first order polynomial that preserves the cell average of v j,h (x). Here, the sub-index p denotes the projected polynomial using (35). We remark that we have used again v h to denote the fluctuation with respect to the stationary solution as in (29). Likewise, x j,c denotes the centre of the element Ω j and the time dependency has been dropped to simplify the notation. This projected polynomial is then used to rewrite the new polynomial fluctuation on the element Ω i as follows: where The normalized weights for the WENO convex combination in (36) follows the classical WENO procedure: The parameters γ l used in this paper are γ 1 = 0.998, and γ 0,2 = 0.001; is set to = 10 −6 . The employed values for r are in {0.25, 2}, and will be better detailed in the next Subsection. Finally, the smooth indicators β l,i are defined as Again, we stress that the general limiting procedure described here is applied once per time-step for either Runge-Kutta or ADER time discretisation.

Double-Loop Limiting Procedure
Even though the component-wise WENO technique is essentially non-oscillatory, the resulting scheme may still generate spurious oscillations, especially at very high order. We, therefore, propose the following strategy.
• First, the procedure (35)- (38) with the value r = 2 is used to produce a candidate solution u c h (x). This candidate solution is evaluated once more by the detector described in Section 4.1.
• Then, if the element is once again considered as troubled, the limiting strategy (35)- (38) will be repeated changing the projector operator. Now, the fluctuation is projected into a constant value given by Additionally, the parameter r will be set to r = 0.25. This strategy will ensure a smoother discrete solution for highly oscillating regions near strong discontinuities.
That finishes the limiting technique description.

Remark 4.
It is essential to notice that the stationary solutions considered in this work are smooth, and therefore the detector does not mark cells as troubled. However, let us suppose that roundoff errors or local extremes provoke that an element is marked as troubled in the presence of a stationary solution. In that case, the use of the fluctuation (29) will ensure that the stationary solution is preserved, since the procedure (35)- (38) is exact for the null operator.

Numerical Test
In this Section, some numerical experiments are conducted to test the proposed numerical first, second, third and fourth-order schemes for Runge-Kutta and ADER time discretisations.
We consider four systems: Burgers, Shallow water, Euler with gravity and the hyperbolic dispersive Gavrilyuk-Favrie equations. For each system, a brief description is given. Then, we test the ability of the numerical scheme to preserve stationary solutions. Afterwards, some numerical tests consisting of a perturbation of a stationary solution are shown alongside the behaviour of the well-balanced limiter technique proposed in this paper.
The CFL number in (11) is always set to 0.9, meanwhile the limiter value M defined in (33) is set to 1, except for the problem (63) that is set to 0.1. For the rest of the limiter parameters appearing in Section 4, we maintain the same values that were defined in that Subsection.
Remark that for the two-spatial domain dispersive system, the CFL condition (11) for DG methods in the case of two dimensional domains reads and we set CFL = 0.9. Besides, for the two dimensional system, the values of the γ weights for the limiter are set to γ c = 0.996 for the central element, and γ = 0.001 for the four neighbouring elements.

Burgers' Equation
Let us consider the one-dimensional scalar balance law consisting of Burgers' equation with a source term written in the form (1) with whose stationary solutions are given by and with H(x) a known continuous function.

Preservation of a Stationary Solution
In this numerical experiment we set as initial condition a stationary solution of the form (43) to ensure the well-balancing properties of the proposed methods for the system (42). The initial condition considered is given by The computational domain is Ω = [−1, 1]. Dirichlet and open boundary conditions were imposed at the left and at the right of the domain, respectively. Figure 1 shows the error between the computed numerical solutions and the stationary one given by e H(x) . As it can be seen in Figure 1 and Table 1, the well-balanced methods can preserve the stationary solution at a final time t = 10 s.

Perturbed Stationary Solutions
We consider now a perturbation of the stationary solution (44) given by The computational domain Ω = [−1, 1] is covered with a coarse mesh of N x = 100 elements. Dirichlet and open boundary conditions are imposed at the left and the right of the domain, respectively. The numerical experiment is run up to time t = 10 s with the fourth-order DG method and for both Runge-Kutta and ADER time discretisations. Figure 2 shows the error between the computed numerical solutions and the stationary one given by e H . As it can be seen, once the perturbation has left the domain, the wellbalanced numerical scheme is able to recover the stationary solution e H , in contrast with the non well-balanced scheme.

Limited Simulation
In this final example, we aim to check the robustness of the proposed limiter technique described in this paper. To do that, we consider the initial condition The computational domain Ω = [0, 2π] is covered with a coarse mesh of N x = 80 elements. Periodic boundary conditions were imposed everywhere. Figure 3 shows the numerical approximation of u with different time-marching discretisations at the final time t = 1.5 s, as well as a reference solution that has been computed by a first order scheme with N x = 20,000 elements. As it can be seen, the proposed schemes provide an accurate approximation without producing any spurious oscillations, even for the case of coarse meshes.

Well-Balanced and Non-Well-Balanced Schemes Comparison
We perform now a comparison between the well-balanced and non-well-balanced schemes for a solution which is not a perturbation of an equilibrium state. Particularly, the following initial condition is considered: Periodic boundary conditions and 100 discretization points are considered in a computational domain Ω = [−1, 1]. Figure 4 depicts the results for a Runge-Kutta second and third order scheme. The reference solution has been computed using a first order method with 20,000 discretization points. As it can be seen, the results shows that the well-balanced and non well-balanced methods are in fact very similar regardless of the order of the method. Similar results are obtained for ADER time marching methods.

Shallow Water Equations
Let us consider the one-dimensional shallow water equations written in the form (1) with where h = h(x, t) is the total water depth, u is the depth-averaged velocity on the x−direction, and H = H(x) is the known still water depth that we will suppose to be continuous. Finally, g = 9.81 is the gravitational acceleration. In this case we adapt the wellbalanced strategy described in Section 3 to preserve the stationary solutions corresponding to water at rest: where η = h − H is the free-surface elevation. A well-balance method could also be design following the ideas described in [93].

Order of Accuracy Test
We begin by checking numerically the order of accuracy of the numerical scheme for a Runge-Kutta and ADER time discretisation for a second-, third-, and fourth-order schemes. As pointed before, we consider here only a well-balanced method that preserves the water at rest solution adapting the strategy described in Section 3. The computational domain is Ω = [−5, 5], and periodic boundary conditions are imposed everywhere. We then simulate the propagation of a Gaussian wave, with the initial condition given by and until the final time t = 1 second, with a set of refined meshes. Table 2 shows the L 1 errors and the numerical convergence rates. The expected convergence order of the Runge-Kutta and ADER DG methods can be observed.

Preserving Stationary Solutions
In this experiment, we set as initial condition a stationary solution of the form (49) to ensure that the well-balanced methods are able to preserve it up to machine precision, while the non well-balanced methods fail. The initial condition considered is given by The computational domain is Ω = [−1, 1] and Dirichlet boundary conditions are imposed everywhere. As it can be seen in Figure 5 and Tables 3 and 4, the well-balanced methods are able to preserve the stationary solution at a final time t = 10 s for both timemarching discretisations and relative coarse meshes. Again, the numerical solutions are compared with the stationary solution (51) as reference.

Perturbed Stationary Solutions
We consider now a perturbation of the stationary solution (51) given by The computational domain Ω = [−1, 1] is covered with N x = 100 elements. Dirichlet and free-outflow boundary conditions are imposed at the right and the left of the domain, respectively.
As it can be seen in Figure 6, the well-balanced scheme can recover the stationary solution at a final time t = 10 s with the fourth-order DG method and for both Runge-Kutta and ADER time discretisations. In this case, we plot the difference between the numerically computed free surfaces and the exact stationary (51), which is zero. Moreover, in order to show the better performance of a well-balanced method against a non well-balanced method, we propose a similar numerical test on a very coarse mesh composed by 20 elements. Here, we consider the initial condition given by The numerical computation is then run up to the short time t = 0.5 s; the perturbation does not reach boundaries in this case, and therefore periodic boundary conditions were imposed everywhere.
As the perturbation is of the order (∆x) N+1 , N + 1 = 4 being the order of the scheme, we can expect that only the well-balanced method yields satisfactory results. In contrast, the non well-balanced scheme introduces oscillations of that order that destroy the solution.
Here we only show the numerical results in Figure 7 for the Runge-Kutta time discretisation. Similar results can be observed for ADER time discretisations.

Limited Simulation
In order to check the robustness of the limiting procedure described in Section 4.3, we consider now a situation far from the equilibrium, where the use of some limiting technique is mandatory. The computational domain Ω = [−5, 5] is covered now with a very refined mesh of N x = 400 elements. The initial condition is set to Periodic boundary conditions are imposed to obtain more interactions between waves. Some snapshots of the computed free surface are depicted in Figure 8. As it can be observed, there are no oscillations at the neighbourhood of the shock waves, which indicates that the limiting approach provides a high overall quality of the computed numerical solution. The solution is depicted alongside a reference solution that has been computed by a first order scheme with N x = 20,000 elements.
Note that the small oscillations appearing with the ADER method could be prevented with an ad hoc tuning of the parameters associated with the limiter described in Section 4. Here, we have maintained the same parameters for both approaches to compare both Runge-Kutta and ADER methods. . Computed free surface for the problem (54) (shallow water equations) with the fourth-order well-balanced Runge-Kutta (left) and ADER (right) numerical methods at times t = 4, 5.5, 6.5 s (from upper to lower panels). An uniform Cartesian mesh of N x = 400 elements has been used.

Compressible Euler Equations with Gravitational Force
We consider the gas dynamic Euler equations written as in (1): ρ ≥ 0 being the density, u the velocity, E the total energy, and H(x) a given continuous gravitational potential. The internal energy e is given by ρe = E − 1 2 ρu 2 . The pressure p ≥ 0 is given by e through the equation of state (ideal gas), where γ > 1 is the adiabatic constant that will take the value 5 3 unless stated otherwise. In this work, we will focus on the hydrostatic steady states solutions given by According to [35], a family of stationary solutions of (55) with u = 0 are where α is a given continuous function and β is given by The numerical methods are then chosen to preserve one of those families given by the choice of α : that results in the following set of stationary solutions with u = 0,

Preserving Stationary Solutions
In this experiment, we set as initial condition a stationary solution of the form (60) to check the well-balancing properties of the proposed methods for the system (55). The initial condition is given by The computational domain is Ω = [−1, 1] and Dirichlet boundary conditions were imposed everywhere. Figure 9 shows the error between the compute numerical solutions and the stationary one given by (61). As it can be seen in Figure 9 and Table 5, the wellbalanced methods can preserve the stationary solution at a final time t = 10 s. Table 5. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (61) (Euler equations).   methods. An uniform Cartesian mesh of N x = 100 elements has been used. Upper and lower panels stands for Runge-Kutta and ADER time-marching discretisations, respectively.

Limited Simulation
We check again the well behaviour of the proposed mixed-order technique described in this paper. To do that, we set H = 0, γ = 1.4, and consider the following discontinuous initial condition given by Then, we set the values for the two well-known numerical tests: • The sod problem [94]: • The Lax problem [95]: ρ L u L p L = 0.445, 0.698, 3.528 , ρ R , u R , p R = 0.5, 0, 0.571 . (64) In both cases, the computational domain Ω = [−5, 5] is covered with N x = 300 elements. The limiter value M defined in (33) is set to 0.01, and the CFL number to 0.5. Free outflow boundary conditions were imposed everywhere. Figures 10 and 11 show the results with the fourth order Runge-Kutta and ADER DG methods for the Sod and Lax problems respectively. The numerical results are compared against a reference solution that has been computed using a first order method with N x = 20,000 elements.

Dispersive Gavrilyuk-Favrie Equations
Finally, we will consider a hyperbolic dispersive shallow water system that was firstly introduced in [96] for the case of a flat bottom bathymetry, and later extended for the nonflat bottom case in [78]. The system can be seen as a hyperbolic relaxation approach of the elliptic Serre-Green-Naghdi equations, and is formally deduced from a master Lagrangian formalism. The system for the case of two-spatial domains reads where Here, h = h(x, y, t) is the water depth, H = H(x, y) is the known still water depth that we will assume to be continuous, and u, v, w are the depth averaged velocities at the x, y and z direction respectively. Furthermore, p(h, ξ) is the non-hydrostatic pressure 81 is the gravitational acceleration, and λ is a large pressure constant parameter that makes the system (66) formally converges to the elliptic Serre-Green-Naghdi equations [97] when λ → ∞. Finally, ξ(x, y, t) is an auxiliary equilibrium variable introduced in the relaxation procedure. In the following, we will set λ = 300 for the numerical test considered.
The system satisfied an extra energy conservation law, and the eigenstructure is fully well known. For more details, the reader is refereed to [78]. As for Shallow-water equations, we will consider the family of stationary lake at rest solutions given by where η = h − H is the free-surface elevation. There are many crucial reasons for considering the numerical validation of the proposed methods for that system here: we will validate the proposed numerical methodology for dispersive shallow flows numerical solutions, and the case of two spatial dimensions; the system introduces a non-homogeneous geometrical source term G(U); the usual steady-states lake at rest solutions for this system contains nonlinear relations between the conserved variables, and that makes standard well-balancing techniques no longer valid to preserve them. Finally, up to our knowledge, this is the first time that a well-balanced high-order method has been proposed for this particular dispersive hyperbolic system.

Preserving Stationary Solutions
In this first numerical experiment, we set as initial condition a stationary solution of the form (67) to check the well-balancing abilities of the proposed methods for the system (66). The initial condition considered (see Figure 12) is given by The computational domain is Ω = [−5, 5] × [−5, 5] and periodic boundary conditions are imposed. As it can be seen in Figures 13 and 14 and Table 6, the well-balanced methods are able to preserve the stationary solution at a final time t = 10 s for both time-marching discretization and relative coarse meshes. Note that the solutions are compared against the stationary solution (68) as reference.  Table 6. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (68) (dispersive Gavrilyuk-Favrie equations).   √ u 2 + v 2 + w 2 for the stationary problem (68) (dispersive Gavrilyuk-Favrie equations) at time t = 10 s with the fourth order non well-balanced (left) and well-balanced (right) methods. An uniform Cartesian mesh of N x × N y = 50 × 50 elements has been used. Upper and lower panels stands for Runge-Kutta and ADER time-marching discretisations, respectively.

Limited Simulation
Next, we check the behaviour of the proposed mixed-order technique described in this paper. To do that, we set the following circular dam-break initial condition: u = v = w = 0, hξ = h 2 , H = 1 − 0.2e −2x 2 , h = H, for x 2 + y 2 ≥ 1 0.25 + H, for x 2 + y 2 < 1.
The computational domain Ω = [−5, 5] × [−5, 5] is covered with a mesh of N x × N y = 200 × 200 elements. We set periodic boundary conditions everywhere. Some snapshots of a section at y = 0 of the computed free-surface are depicted in Figure 15 with a fourthorder Runge-Kutta time marching discretization. As it can be observed, the comparison between a limited and a non-limited simulation shows that there is no oscillations at the neighbourhood of the dispersive shock waves. Note that the high-oscillatory and overamplified waves appearing in Figure 15 from the unlimited numerical solution, are not of the nature of this hyperbolic dispersive system.

Conclusions
We have presented a novel and general technique for preserving stationary solutions for Discontinuous Galerkin methods. This general approach is able to transform standard DG schemes in exactly well-balanced arbitrary high-order Discontinuous Galerkin methods for both Runge-Kutta and ADER time-marching discretisations.
For continuous source terms, the methodology only uses a well-suited definition of a well-balanced reconstruction operator and not on the Riemann-solver. Therefore, the well-balanced technique does not rely on the particular details of the system. That makes the method very general and can be applied to any system of hyperbolic PDEs, as long as the desired stationary solutions to be maintained are known.
Furthermore, the standard WENO limiting procedure is carefully considered. Here, we extend classical ideas to propose a novel WENO limiter that preserve stationary solutions. That makes the final numerical methods proposed here arbitrary high-order, essentially non-oscillatory and well-balanced.
We have tested the numerical methodology derived here with some standard systems, proving the well-balancing properties of the approach, the high-order accuracy, and the general well behaviour of the technique in the presence of strong discontinuities.