On High Order ADER Discontinuous Galerkin Schemes for First Order Hyperbolic Reformulations of Nonlinear Dispersive Systems

This paper is on arbitrary high order fully discrete one-step ADER discontinuous Galerkin schemes with subcell finite volume limiters applied to a new class of first order hyperbolic reformulations of nonlinear dispersive systems based on an extended Lagrangian approach introduced by Dhaouadi et al. (Stud Appl Math 207:1–20, 2018), Favrie and Gavrilyuk (Nonlinearity 30:2718–2736, 2017). We consider the hyperbolic reformulations of two different nonlinear dispersive systems, namely the Serre–Green–Naghdi model of dispersive water waves and the defocusing nonlinear Schrödinger equation. The first order hyperbolic reformulation of the Schrödinger equation is endowed with a curl involution constraint that needs to be properly accounted for in multiple space dimensions. We show that the original model proposed in Dhaouadi et al. (2018) is only weakly hyperbolic in the multi-dimensional case and that strong hyperbolicity can be restored at the aid of a novel thermodynamically compatible GLM curl cleaning approach that accounts for the curl involution constraint in the PDE system. We show one and two-dimensional numerical results applied to both systems and compare them with available exact, numerical and experimental reference solutions whenever possible.


Introduction
Nonlinear dispersive systems can be found in many different areas of computational mechanics, ranging from large scale dispersive free surface shallow water flows [12,28,62,87,93,94] over multi-phase flows with surface tension [13,39] down to quantum fluid mechanics [14,15,65,75]. A common difficulty in the above applications is that the nonlinear timedependent governing partial differential equations (PDE) typically contain either higher order spatial and temporal derivatives, or a subset of elliptic equations, and therefore the integration with simple explicit finite volume and discontinuous Galerkin finite element schemes becomes unfeasible.
Therefore, very recently, there has been an increasing interest in rewriting nonlinear PDE systems with higher order derivatives under the form of nonlinear hyperbolic relaxation systems, which contain at most first order derivatives in space and time, in conjunction with potentially stiff algebraic relaxation source terms. The idea goes back to the seminal work of Cattaneo [25], who proposed a hyperbolic reformulation of the parabolic heat equation. More recent work on the topic also regards the hyperbolic reformulation of advection-diffusion equations [80,84,85,99], of the compressible Navier-Stokes equations [17,47,88] as well as hyperbolic reformulations of nonlinear dispersive systems [3,4,26,51,52,63,64,79,90]. We would like to point out a particularity of the hyperbolic dispersive system proposed in [3], since it can be derived from the depth averaged compressible Euler equations, while most depth-averaged shallow water systems are usually derived from the governing equations of an incompressible fluid. Also the well-known theory of rational extended thermodynamics [82] makes use of hyperbolic relaxation systems to model the effects of higher order derivative terms.
In this paper we focus in particular on the hyperbolic reformulation of the Serre-Green-Naghdi model introduced by Favrie and Gavrilyuk in [56], as well as on the hyperbolic reformulation of the nonlinear defocusing Schrödinger equation proposed by Dhaouadi et al. in [38]. Both systems were rigorously derived from an extended Lagrangian formalism and therefore have a rather similar mathematical structure. The hyperbolic reformulation of the nonlinear Schrödinger equation [38] has the additional difficulties that it is only weakly hyperbolic in multiple space dimensions and that it is endowed with a curl involution constraint that needs to be properly accounted for. In this paper we will make use of the hyperbolic generalized Lagrangian multiplier (GLM) curl cleaning approach of [27,46] that goes back to the hyperbolic GLM divergence cleaning technique of Munz et al. for Maxwell and MHD equations, see [35,83]. In alternative to the GLM method proposed here, also exactly curl-preserving schemes could be used, see e.g. [1,11,27,66,67,71], but they require an appropriately staggered mesh and are therefore not as easy to implement in an existing general purpose DG solver as the simple GLM method, which only requires the solution of additional PDEs for the cleaning quantities. The main novelty proposed in the present paper is a new GLM curl cleaning that is also thermodynamically compatible with the conservation of total energy.
To integrate the governing PDE systems under consideration in this paper we will make use of high order accurate explicit discontinuous Galerkin (DG) finite element schemes. The DG framework for hyperbolic conservation laws goes back to the seminal work of Cockburn and Shu [29][30][31][32] and was later generalized to hyperbolic equations with parabolic terms in [5,6,33,34]. The first extension of DG schemes to dispersive PDE of the Kortevegde-Vries (KdV) type with up to third order spatial derivatives was achieved in [106,107], while DG schemes for PDE with even higher order spatial derivatives were first tackled in [70]. Discontinuous Galerkin finite element schemes for the solution of nonlinear dispersive Boussinesq-type equations were forwarded in [50,54,55]. To overcome the very severe time step restriction of explicit DG schemes applied to dispersive equations, which require the time step to scale with the cube of the mesh spacing (Δt ∝ Δx 3 ), see [106], fully implicit spacetime DG schemes for dispersive problems were introduced in [45]. However, although the resulting schemes are unconditionally stable they are computationally very expensive due to the ill-conditioned algebraic systems that must be solved in each time step. For this reason, in this work we prefer the numerical solution of hyperbolic reformulations of nonlinear dispersive systems, which allow the straightforward use of simple explicit DG schemes for hyperbolic PDE with the usual CFL-type time step restriction, according to which the time step must be chosen proportional to the mesh spacing (Δt ∝ Δx), rather than more complex implicit DG schemes for PDE with higher order derivatives. In alternative to high order DG methods, also high order residual distribution (RD) schemes can be applied to Boussinesqtype equations, see [90].
The rest of this paper is structured as follows: in Sect. 2 we present the two nonlinear dispersive models that we want to study, namely the hyperbolic reformulation of the Serre-Green-Naghdi system of dispersive free surface water waves forwarded by Favrie and Gavrilyuk in [56] and the hyperbolic reformulation of the nonlinear defocusing Schrödinger equation proposed by Dhaouadi et al. in [38]. Compared to the original work on the hyperbolic reformulation of the Serre-Green-Naghdi system [56] in this paper we also add the bottom slope term and design an exactly well-balanced numerical scheme for the resulting system that is capable of preserving stationary lake at rest solutions exactly at the discrete level for arbitrary bottom topography. Compared to the work on the hyperbolic reformulation of the Schrödinger equation [38] in this paper an additional term is added to the system in order to restore Galilean invariance and also the curl involution constraint of the system is explicitly taken into account, which is necessary for the multi-dimensional case. As already mentioned previously, in this paper we will make use of a hyperbolic generalized Lagrangian multiplier (GLM) approach, see [27,35,46,83], which also restores strong hyperbolicity of the system. In Sect. 3 we present the framework of arbitrary high order derivative (ADER) discontinuous Galerkin finite element schemes used for the numerical solution of the hyperbolic governing PDE systems. Section 4 is devoted to the presentation and discussion of the numerical results obtained for both systems. Wherever possible, we compare with exact, numerical or experimental reference solutions. In Sect. 5 we draw some conclusions and give an outlook to future research.

Governing Equations
Throughout this paper we will make use of the Einstein summation convention over two repeated indices. We further denote the time with t and the Cartesian coordinate axes by x k for k ∈ {1, ...d} with d the number of space dimensions. Sometimes we will also make use of the notation x := x 1 , y := x 2 and ∂ x = ∂ ∂ x .

A Hyperbolic Reformulation of the Serre-Green-Naghdi Model
The governing PDE system generalizing that proposed in [56] to reformulate the Serre-Green-Naghdi model as a first order hyperbolic system and derived from the extended Lagrangian variational principle reads as follows (see Appendix A for details): with the nonhydrostatic pressure , the gravity constant g = 9.81 m 2 /s and where the parameter λ having the dimension m 2 /s 2 is a large (compared to the squared velocity of surface waves gh) free parameter which makes the system (1) -(5) tend to the original Serre-Green-Naghdi model in the limit λ → ∞. The rigorous proof of this fact has been established in [40] in the case of flat bottom. Note that for λ → ∞ the term (1 − η/h) → 0 since η → h, hence the product will remain finite. Compared to [56] in system (1)-(5) the mild bottom slope terms are added (see Appendix A for details).
The system admits the energy conservation law which is a direct consequence of the Noether theorem.
For a hyperbolic reformulation of the Serre-Green-Naghdi model without mild bottom assumption and directly derived from the compressible Euler equations, the reader is referred to [3]. The derivation from the compressible Euler equations presented in [3] can also be used to establish a connection between the square of the sound speed in the compressible medium and the penalty parameter λ used in the model (1)- (5). However, the model proposed in [3] is non-conservative even for flat bottom (see also [51,52]). In the rest of this paper we sometimes refer to the model (1)-(5) also as the Favrie-Gavrilyuk (FG) model, since it is the extension of [56] to the case of mildly varying bottom.

Eigenstructure
The hyperbolicity of system (1)-(5) was already studied in [56] for the one dimensional case with flat bottom, while here we show the results for the two-dimensional situation with variable bottom. The eigenvalues of system (1)-(5) in x 1 direction are with the celerity a = gh + 1 3 λ η 2 h 2 and the associated right eigenvectors and with the auxiliary quantities

Linear Dispersion Analysis and Comparison with Other Hyperbolic Models
Following [3,51,73,77], we carry out a linear dispersion analysis of the system (1)-(5) for flat bottom, ignoring all terms involving b. To simplify the expressions, we consider only the one-dimensional case and introduce the vector of primitive variables V = (h, u, η, w) T . Rewriting of the governing PDE system in terms of the primitive variables leads to with After decomposition of the primitive variables into a stationary component V 0 plus a small time-dependent fluctuation V and assuming that S(V 0 ) = 0, we get the linearised system that can be studied using standard Fourier analysis. The ansatz V (x, t) =Ve i(κ x−ωt) with i 2 = −1 leads to Hence, from (11), we obtain In particular, taking the rest state V 0 = (H , 0, H , 0) the linearised matrix related to system (1)-(5) in primitive variables reads The eigenvalue problem M 0 V = ωV can be solved numerically, thus obtaining the phase velocity C p = ω κ . For comparison with other available models in the bibliography, it is important to note that the parameter λ is chosen as λ = 3c 2 = 3α 2 g H with c an artificial sound speed, see [3,51]. The relative error against the phase velocity given by the linear theory of Stokes, , is depicted in Fig. 1 for κ H ∈ [0, 3]. The errors obtained for the system in [56] given by (1)-(5) and the hyperbolic Serre-Green-Naghdi model for general bottom topographies proposed in [3] match pretty well, as expected, since both systems are hyperbolic reformulations of the Serre-Green-Naghdi model. The relative errors for the Madsen and Sørensen model [76,77], which has improved dispersion characteristics, are also included for comparison. Moreover, the right plot in Fig. 1, shows the sensitiveness of the model to variations of λ. We observe that the error decreases when increasing the value of the parameter and the error curves almost overlap for large enough values of λ.

Hyperbolic Reformulation of the Nonlinear Schrödinger Equation
The first order hyperbolic reformulation of the nonlinear Schrödinger equation according to [38] reads with the total pressure P = 1 2 ρ 2 − 1 4ρ p m p m + λ η 1 − η ρ and where λ 1 and β 1 are two free parameters which make the system (15)- (19) tend to the original hydrodynamic model of the Schrödinger equation in the limit λ → ∞ and β → 0. Please note that compared to [38] the governing equation for p is different in this paper, since the term v m has been added in order to restore the Galilean invariance of the system, see [27,47,60,92]. We note that the mathematical structure of system (15) is very similar to the one of (1)-(5), apart from the additional field p k that was not present in the hyperbolic reformulation of the Serre-Green-Naghdi model. We stress that system (15) is endowed with the curl involution constraint ∇ × p = 0, which is a direct consequence of (15) if the constraint is satisfied by the initial data. Recall that with p = ( p 1 , p 2 , p 3 ) T we denote the vector with components p k . In [38] the hyperbolicity of (15) was only studied in the one-dimensional case. In this paper we present a complete study in three space dimensions, which yields some interesting findings.

Eigenstructure of the Original Weakly Hyperbolic System
The eigenvalues of system (15)- (19) in x 1 direction are (20) with the celerities c λ = ρ + 1 The associated right eigenvectors read with the auxiliary quantities As one can easily see, two eigenvectors are missing, hence the system (15)-(19) based on [38] is only weakly hyperbolic in the multi-dimensional case. The reason behind this is that the curl involution constraint has not yet been taken into account. In order to fix this issue, in Sect. 2.2.3 we will introduce an augmented system that accounts for the curl involution on p via a hyperbolic generalized Lagrangian multiplied (GLM) curl cleaning approach, which goes back to the ideas forwarded by Munz et al. [35,83] for divergence-type constraints and in [27,46] for curl-type constraints.

Stationary Equilibrium Solution in Cylindrical Coordinates
In this section, we derive a two-dimensional stationary equilibrium solution of system (15)- (19). For this purpose, the equations are first rewritten more conveniently in cylindrical coordinates r − φ − z, with the velocity vector written as v = (v r , v φ , v z ) and the vector field p = ( p r , p φ , p z ). A substantially simplified system in radial direction is obtained by assuming ∂ φ = ∂ z = 0 and v z = p z = 0. The final system in radial direction reads as follows: where in the last two equations the constraint ∇ × p = 0 has been used. We are now looking for a stationary solution, hence ∂ t = 0. We furthermore assume a vortex-type solution with v r = 0 and p φ = 0 so that v · p = 0. With these hypothesis from (25) we obtain w = 0. The above system (22)-(28) thus reduces to the following ODE system in radial direction: We now prescribe a radial profile for the density ρ and for the radial component p r as follows: and From (30) one then directly obtains the profile η(r ) as while the angular velocity profile v φ (r ) is obtained from (29) and reads The exact expressions for the profiles of η(r ) and v φ (r ) can be easily obtained via a computer algebra system and are thus not explicitly reported here.

Augmented GLM Curl Cleaning System
Following the original ideas of Munz et al. [35,83] concerning hyperbolic divergence cleaning for the Maxwell and MHD equations and their extension to curl involutions in [27,46], we propose the following augmented system, where eventual curl errors produced by the numerical scheme are transported away by a Maxwell-type subsystem: with the curl cleaning speed a c = const > 0. The additional GLM curl cleaning terms have been highlighted in blue, for convenience. From the last two equations we recover the involution constraint ∇ × p → 0 in the limit a c → ∞. The equations are also compatible with the energy conservation law (see Appendix B), hence the hyperbolic GLM cleaning proposed in this paper is a thermodynamically compatible one, in contrast to the previous hyperbolic GLM cleaning approaches proposed in [27,35,37,46] that were not compatible with the conservation of total energy. Note, however, that the GLM divergence cleaning for MHD proposed in [37] was compatible with the conservation of mathematical entropy.

Eigenstructure of the Augmented GLM System
The eigenvalues of the thermodynamically compatible augmented GLM curl cleaning system (35)-(40) in x 1 direction are The associated right eigenvectors are with the auxiliary quantities For p 2 = 0 or p 3 = 0 the eigenvectors r 3,10 are linearly independent of the others and thus the system is strongly hyperbolic. It becomes weakly hyperbolic only for p 2 = p 3 = 0. One immediately realizes that the eigenvectors r 1,12 , r 2,11 and r 5,6,7 coincide exactly with those of the original system without GLM curl cleaning.

Numerical Method
Both dispersive models (1)- (5) and (35)- (40) can be cast into the following general form with the state vector Q ∈ Ω Q ⊂ R m , the state space Ω Q , the flux tensor For the numerical simulation of the hyperbolic reformulations of the nonlinear dispersive systems introduced in the previous section, and which can all be cast into the general form (44), we employ the family of high order accurate fully-discrete one-step ADER discontinuous Galerkin schemes with a posteriori subcell finite volume limiter, see e.g. [17,42,43,49,108]. In the following we provide a brief description of the method. For further details the reader is referred to the above references.

Fully Discrete One-Step ADER-DG Schemes
The governing PDE systems (1)- (5) and (35)- (40) are discretized on a computational domain Ω employing a uniform Cartesian grid with elements x and y direction, respectively. Denoting by u h (x, t n ) the discrete solution of (44) written in the space of piecewise polynomials of degree N , the discrete solution within Ω i is sought under the form Here, φ l (x) = ϕ l 1 (ξ )ϕ l 2 (η) are the basis functions, which are tensor products of onedimensional basis functions ϕ l m (χ) on the unit interval χ ∈ Ω ref = [0, 1]. The mapping from the reference coordinates 0 ≤ ξ, η ≤ 1 to the physical coordinates reads x = x i − 1 2 Δx +ξΔx and y = y i − 1 2 Δy + ηΔy. Here, l is a multi-index, referring to the one-dimensional basis functions ϕ l m to be used in the tensor product. We employ Lagrange interpolation polynomials passing through the Gauss-Legendre quadrature points of a Gaussian quadrature formula with N + 1 nodes. This leads by construction to an orthogonal basis. Due to the nodal tensor-product basis, the scheme can be written in a dimension-by-dimension fashion and all operators can be decomposed into products of one-dimensional operators. Multiplication of the governing PDE system (44) by test functions φ k , which are identical to the basis functions, and integration over a space-time control volume Using (45) and integrating the flux divergence term by parts in space and the time derivative by parts in time the above weak problem becomes ⎛ with the maximum wavespeed at the interface s max = max(|λ k (q − h )|, |λ k (q + h )|) and the matrix I wb > 0. For the hyperbolic reformulation of the Schrödinger equation I wb = I is simply chosen equal to the identity matrix, while for the hyperbolic Serre-Green-Naghdi model (1)-(5) I wb needs to be chosen in a special manner in order to make the numerical scheme well-balanced for the lake at rest solution [10,59,69], see Sect. 3.3. In alternative, also the more sophisticated generalized Osher-type scheme forwarded in [48] can be used instead of (48). Here, q − h and q + h denote the boundary-extrapolated values of the space-time predictor from within the element and its neighbor, respectively. The jump terms in the non conservative products at the boundaries are treated via a path conservative scheme as forwarded by Castro, Parés and collaborators in [21,22,24,81,86], which are based on the theory of Dal Maso, Le Floch and Murat [78] on nonconservative hyperbolic PDE systems. For a more detailed discussion on the topic, see also [23] and references therein. Path-conservative schemes were generalized to higher order DG schemes for the first time in [43,89]. Based on the pathconservative framework, the construction of so-called well-balanced schemes for shallow water type models is rather straightforward. The jump term in the non-conservative product is computed via a path integral between the two boundary extrapolated values related to the face, q − h and q + h as where we use the simple straight line segment path In this paper, the path integral (49) is approximated via a simple trapezoidal quadrature rule, which is enough to obtain a well-balanced scheme, as we will see later in Sect. 3.3. In order to achieve an essentially non-oscillatory behaviour at discontinuities and in the regions of strong gradients, the ADER-DG schemes used in this paper are supplemented with an a posteriori subcell finite volume limiter, as detailed in [49,108].

Local Space-Time Predictor
The local space-time predictor q h (x, t) is obtained via a weak formulation of the governing PDE system in space-time, as proposed in [42,44]. This allows to avoid the cumbersome Cauchy-Kovalewskaya procedure used in the original ADER finite volume schemes of Toro et al. [20,96,97,100,101]. The element-local space-time predictor solution is introduced as follows: with the space-time basis functions θ l = θ l (x, t) = ϕ l 0 (τ )ϕ l 1 (ξ )ϕ l 2 (η), which are tensor products of the spatial basis functions already introduced before and an additional nodal basis function for the time dependency, and with t = t n + τ Δt. Multiplying (44) by a space-time test function θ k and integrating over Ω i × t n , t n+1 , one obtains Integration by parts of the first term in time yields (53) which is an element-local system for the unknown degrees of freedomq k of the spacetime predictor q h (x, t) and which can be computed in terms of the known spatial degrees of freedomû n l,i of the discrete solution u h (x, t n ). The solution of (53) can be found via a fast-converging iterative fixed point scheme, the convergence of which was proven in [17]. As any explicit numerical method for hyperbolic systems, our scheme is subject to a usual CFL-type time step restriction, requiring with C < 1, see [42], and where the maximum is taken over all cells and degrees of freedom in the computational domain. It is clear that (54) contains the parameter λ, but the time step scales only linearly in Δx and not cubically. Furthermore, the proposed explicit DG scheme does not require the solution of big linear systems, unlike the fully implicit space-time DG scheme for dispersive systems proposed in [45].

Well-Balanced Property for the Favrie-Gavrilyuk Model with Variable Bottom Topography
In this section we prove the well-balanced property [69,104] of the first order version of our ADER-DG schemes in one space-dimension, i.e. for N = 0 and d = 1, which remains the key ingredient of the well-balancing also for higher order schemes since it involves the numerical flux and the jump terms. For high order well-balanced reconstructions see [21,22,57,58]. Later in this section, we also study the well-balanced behavior of the actual implementation of the numerical scheme. For N = 0 we have only one single test and basis function φ 1 = θ 1 = 1 and one single degree of freedom in space and timeq n 1,i =û n 1,i = Q n i per cell, which corresponds to the usual cell average Q n i . The scheme (47) then reduces to the following path-conservative finite volume method: with the numerical flux and the path-conservative jump term that in one space dimension and For the system (1)-(5) the following well-balanced identity matrix has been chosen throughout this paper: Lake at rest solutions for system (1)-(5) are characterized by and are stationary solutions of (1)-(5) for arbitrary bottom b. We now prove that the scheme (55) maintains solutions of the type (60) exactly at the discrete level for all times. Since η = h and w = 0, it is obvious that the discrete source term S Q n i vanishes. We now analyze the numerical flux G Q n i , Q n i+1 and the path-conservative jump term D Q n i , Q n i+1 . For ζ = const. and η = h we have p = 0 and the following relations for the jumps in h, η and b: From these relations and η = h we can compute the jump in the conservative variable hη as We thus obtain for the jump in the state vector and with p − = p + = 0 and v − 1 = v + 1 = 0 the path integral of B 1 becomes  From (65) and (64) it follows that hence the path-conservative jump term vanishes. Next, we analyze the numerical flux (56). Since The last term to analyze is the numerical viscosity, for which we obtain sinceh = h − + h + = 2h − − Δb and thus −2h − Δb + Δb 2 +hΔb = 0. Therefore, also the numerical flux G q − h , q + h = 0 vanishes at each element interface. Hence, for lake at rest solutions (60) the scheme (55) reduces to i.e. lake at rest solutions are exactly preserved at the discrete level, or, in other words, the scheme is well-balanced.
To verify the well-balanced property also numerically, we carry out the following simple test. In a 1D domain Ω = [−5, +5] we define the following initial condition, with H (x) being the Heaviside function: which is a lake at rest solution (60) of system (1)- (5). The domain Ω is discretized by 200 equidistant control volumes and simulations are run with the first order scheme (55) using (56) with (59) and (57) until a final time of t = 2 using a CFL number of CFL= 0.5 and setting λ = 1200. We report the L ∞ errors for the variables h, v 1 , η and w at time t = 2 in Table 1 using single, double and quadruple precision. The computational results clearly show that the method is well-balanced up to machine precision, as expected.

Numerical Convergence Study
In this section, we simulate where V = √ g(H + A) is the velocity of the soliton and the similarity variable is denoted by ξ . Under these assumptions we get ∂ t Q = −V Q and ∂ x Q = Q . Hence, the PDE system (44) can be rewritten as with the matrix A(Q) = ∂f/∂Q + B 1 (Q) containing the flux Jacobian and the nonconservative product. For flat bottom, as considered here, the latter does not contribute.
The PDE system (44) finally reduces to the following nonlinear ODE system with initial condition Q(ξ 0 ) = (H , 0, 0, H − ε, 0, 0) and with I the identity matrix. For our numerical tests we set g = 9.81, H = 1, A = 0.1 and the perturbation ε = 10 −10 . The ODE system (73) is solved with a 6 th order DG scheme in time, see [41], in order to provide a highly accurate initial condition for the solitary wave of the hyperbolic SGN system (1)- (5). We now simulate the propagation of the solitary wave until t = 2 using ADER-DG schemes of polynomial approximation degrees N = 2, 3, 4, 5 on a sequence of successively refined meshes. The L 2 errors and corresponding numerical convergence rates are shown in Table 2. Overall we find the expected convergence order N + 1 of our high order ADER-DG schemes.
In order to investigate the dependence of the solutions of the model (1)-(5) on the parameter λ more thoroughly, we use the high order ODE solver [41] in order to calculate the shape of a solitary wave moving with velocity V = 4 in still water of depth H 0 = 1 as a function of λ ∈ [75, 6000], i.e. with λ varying over about two orders of magnitude. The dependence of the shape of the resulting solitary waves on λ is depicted in Fig. 2. One can observe that the shapes for λ = 1200 and λ = 6000 match very well. For this reason, if not stated otherwise, in the following numerical experiments we will set λ = 1200, since λ = 6000 does not give visible improvements for practical simulations, but due to the CFL condition reduces the time step by more than a factor of 2 compared to the choice λ = 1200.
To give an idea of the computational cost of the approach presented in this paper, we compare it with the cost of a fully implicit space-time DG scheme [45] applied to the dispersive system of Madsen and Sørensen [76] with third order derivatives. The computational setup is the following: we consider the 1D domain Ω = [−20, +20] with periodic boundary Table 2 Numerical convergence results for high order ADER-DG schemes of polynomial approximation degree N ∈ {2, 3, 4, 5} for the hyperbolic reformulation of the Serre-Green-Naghdi model (1)-(5) according to [56]. A uniform Cartesian mesh composed of N x × 4 elements has been used. The L 2 error norms and corresponding orders os convergence (in bold) refer to the variables v 1 , η and w at a final time of t = 2 conditions in which a solitary wave with H 0 = 1 and V = 4 is moving, see also Fig. 2 for a graphical representation of the shape of the wave. At time t = 10 the exact solution of the problem is given again by the initial condition. We apply the explicit ADER-DG scheme to the hyperbolic model (1)-(5) with λ = 1200 and the implicit space-time DG scheme [45] to the Madsen and Sørensen model [76], which contains higher order derivatives. In both cases solitary waves with the same parameters are considered (H 0 = 1, V = 4). Both schemes are nominally fourth order accurate, i.e. with polynomial approximation degree N = 3 and the computational mesh used in both simulations is composed of N x = 200 elements. In Table 3 we report the obtained L 2 error norms for h and hv 1 at time t = 10 for both schemes, as well as the necessary CPU time and the time per degree of freedom update (TDU). The TDU metric is computed as TDU = WCT/(N t N x (N + 1)), where WCT is the total wall clock time needed for the simulation, N x is the number of elements in space, N + 1 is the number of degrees of freedom per element and N t is the total number of timesteps needed to reach the final time. Both simulations were run in serial on one single CPU core of an Intel i9-7900X workstation with 32 GB of RAM. From the results shown in Table 3 one can see that both methods reach comparable errors in the variable hv 1 , but the explicit ADER-DG scheme applied to the hyperbolic model (1)- (5) is about three times faster compared to the fully implicit space-time DG scheme applied to the dispersive model [76]. The cost per single degree of freedom update (TDU) is three orders of magnitude larger for the fully-implicit  Table 3 Comparison of errors and CPU times for the propagation of a solitary wave with H 0 = 1 and V = 4 at time t = 10 using the explicit ADER-DG scheme applied to the hyperbolic system (1)-(5) with the fully implicit space-time DG scheme described in [45] applied to the Madsen and Sørensen model [76]. We also report the time needed to update a single degree of freedom for one time step (TDU) space-time DG scheme. We expect that in multiple space dimensions and in the context of a parallel implementation on distributed memory machines, the comparison may be even more in favor of the explicit scheme applied to the hyperbolic reformulation of the dispersive system.

Solitary Wave
As second part of the previous test, we now analyse the propagation of a solitary wave of the original Serre-Green-Naghdi system over a flat bottom for a longer simulation time t = 100/ √ g (H + A). The exact solution of a solitary wave for the original Serre-Green-Naghdi equation reads with ξ = x − V t and V = √ g(H + A), see e.g. [56]. As suggested in [38], the remaining variables of the system are set as i.e. with η = h the system is initialized so that p = 0. For this test we set H = 1, A = 0.2 and λ = 1200. The domain is again Ω = [−50, +50] × [0, +25] and periodic boundary conditions are defined everywhere. The results obtained with an ADER-DG scheme of polynomial approximation degree N = 4 on a uniform Cartesian mesh composed of 1080 × 4 elements after one complete revolution of the soliton are depicted in Fig. 3 and are compared against the exact solution given by (74)- (76). As expected, the shape of the soliton has been well preserved, showing only some small spurious oscillations away from the soliton.

Solitary Wave Over a Step
The solitary wave over a step benchmark proposed in [93] and also used in [3,4] is employed as a first test with nontrivial bottom topography aiming to compare the numerical results obtained The still water depth is H = 0.2 and λ = 1200 is chosen. A soliton of amplitude A = 0.0365 is initially located at x sol = −3. To generate the initial conditions for the different variables the corresponding ODE system (73) is solved using a 10 th order DG scheme in time. Periodic boundary conditions are defined in y direction while on the left and right boundaries we set the initial data as Dirichlet boundary conditions. The domain is discretized with 2000 elements in the horizontal direction and the simulation is run with polynomial degree N = 3. A 1D cut of the free surface obtained at times t ∈ {0, 2.148, 4.296, 6.444, 8.592, 10.74} is depicted in Fig. 4. As expected, the soliton initially propagates smoothly growing when approaching the bottom step. Then, it splits into two transmitted waves in forward direction. A secondary reflecting wave, followed by a train of small dispersive waves, is generated and propagated in the opposite direction to the incoming soliton. The use of the a posteriori FV limiter results in the damping of the small spurious oscillations appearing in correspondence to the obstacle that could be seen for the hyperbolic dispersive models used in [3,4]. Figure 5 reports the time evolution of the ratio between the wave amplitude and the still water depth, A H = h+b−H H , at the fixed observation points x ∈ {−9, −3, 0, 3, 6, 9}. The numerical solution of the model (1)-(5) obtained in this paper agrees well with the results obtained using the alternative hyperbolic reformulations of the SGN model employed in [3,4,51]. Moreover, we also observe a good agreement with experimental data.

Periodic Waves Over a Submerged Bar
We now study the results obtained for the periodic waves over a submerged bar benchmark presented in [7,8,72]. The computational domain is depicted in Fig. 6. The physical domain is enlarged in x-direction to introduce damping zones of length L rel = 10 allowing to define a wavemaker and an absorbing boundary condition guaranteeing a smooth transition between the target solution u * and the solution u h inside the computational domain. To define the periodic signal at the left boundary a sinusoidal function is employed, where the amplitude, A ∞ , is selected to match the initial amplitude of the periodic waves given by the experimental data at the wave gauge S 0 = 0 and ω = 2π/T is the angular frequency, T the wave period, also fixed attending to the experimental data. Regarding the absorbing boundary condition we define, as target solution, preventing wave reflection. Inside the relaxation zones the solution in an element Ω i is computed as with d i the distance between the barycentre of Ω i and the corresponding boundary, x L or x R . Periodic boundary conditions are defined in y direction. The experimental data available corresponds to two different settings. The first of them, a low frequency test (LF), is characterized by A ∞ = 0.01, T = 2.0205. Meanwhile, for the high frequency case (HF) we have A ∞ = 1.286E − 2, T = 1.2624. The polynomial degree in both simulations is N = 3. The a posteriori subcell finite volume limiter is employed to prevent spurious oscillations, see [17,49,108]. The parameter λ is chosen as λ = 1200.
To analyse the results obtained, we consider six wave gauges located in S 1 = 10.8,  is available. The numerical results for the wave amplitude, A = h + b − ζ 0 , are reported in Fig. 7 for the low frequency case and in Fig. 8 for the high frequency test. For the low frequency case, two different mesh resolutions and two different values of λ have been used, in order to show that the obtained numerical results have already reached a good level of independence in terms of mesh spacing and λ. For the high frequency testcase, two different values of λ have been used in order to show that the chosen values of λ are large enough for practical simulations. The results are quite satisfactory and comparable with the ones obtained in [3]. The discrepancies, with respect to the experimental data, observed for S 5 and S 6 in LW may be due to the suboptimal linear dispersive properties of the model considered, see Sect. 2.1.2. To improve these results, models with better dispersion characteristics, like the Madsen & Sørensen model [77], may be employed in a future work. Fig. 8 Time evolution of the wave amplitude A(t) for the high frequency periodic waves over a submerged bar test at wave gauges S 1 − S 6 . Black dots correspond to experimental data [8], blue solid lines are the numerical results obtained for model [56] with variable bottom (1)-(5) using an ADER-DG scheme with N = 3, a mesh of 2800 elements and parameter λ = 1200; red dashed lines are the solution for λ = 2400 (Color figure online)

Favre Waves
In this section we consider the numerical simulation of so-called Favre waves, which are dispersive shock waves or undular bores developing when a fluid layer with a free surface hits against a wall, see [51,56] for numerical investigations and [102] for the corresponding experimental studies. The computational domain is Ω = [−180, 0] × [−1, +1] and it is covered by 1200 × 4 ADER-DG elements of polynomial approximation degree N = 3. The initial condition for this test is where F is the relative Froude number of the impact velocity. As suggested in [56], instead of a wall at x = 0 an alternative symmetric impact problem can also be considered. For the subsequent tests we choose the still water depth H = 1 and the parameter λ is set to λ = 1200. In order to capture the effect of wave breaking, only in this section the following simple wave breaking mechanism is employed, see [3,51]. For this purpose we modify the governing PDE for hw by adding a source term as follows: with and the switch function For our numerical tests we set b 1 = 0.5 and b 2 = 0.25. In Fig. 9 we show the computational results obtained for Froude F = 1.15 at a final time of t = 54 using two different meshes, one with 1200 elements in x direction and a finer one with 2400 elements in x direction. We can note an excellent agreement of the two solutions with each other, showing that mesh convergence has been reached for this test. In Fig. 10 we compare our numerical results with the experimental data reported in [56,102]. The upper symbols refer to the amplitude of the first wave, while the lower symbols refer to the amplitude of the trough after the first wave. The agreement between numerical simulations and experimental data is very good up to F = 1.35, which is already in the range of Froude numbers where wave breaking occurs. This means that the simple wave breaking mechanism adopted here is valid at least for moderate Froude number flows with breaking waves.

Solitary Wave Run-Up onto a Plane Beach
Synolakis [95] carried out laboratory experiments for solitary incident waves to study propagation, breaking and run-up over a planar beach with a slope 1 : 19.85. Since then, many researchers used his data to validate numerical models, as in [3,[51][52][53]68,90,91], among others. Accordingly, this test case is here used to assess the capability of the proposed methodology to describe shoreline motions and wave breaking when it occurs. The bathymetry of the problem is described in Fig. 11. As initial condition we use the solitary wave solution of the genuinely SGN system given by (74)- (75) with where p is the non-hydrostatic pressure known for the SGN model. Note that for practical purposes, this analytical solution of the elliptic SGN model can be used as an approximate solution for the hyperbolised system.  [56,102]. For this test the simple wave breaking mechanism (83) has been employed The relaxation speed λ is set to 1200. In order to take into account the glass surface roughness, the usual and simple bottom friction model used for the standard shallow-water equations is adopted via a usual Manning-type friction formula. For this purpose, we modify the governing equation for hv i by adding a source term as follows: The a posteriori subcell finite volume limiter is activated in the presence of spurious oscillations as well as for small values of the total water depth h < 10 −2 . Following [103], we employ a second-order scheme by using a TVD polynomial reconstruction procedure using the MUSCL slope limiter, which takes into account the positivity of the water height, [53], and the well-balancing for equilibrium. The time evolution stage to the half time level is then computed via a standard second-order TVD Runge-Kutta method, [61]. Finally, freeoutflow boundary conditions are considered and the CFL number is set to 0.9. Figure 12 depicts the snapshots at times t ∈ {4.7, 6.3, 7.9, 9.5}. Excellent results are obtained for the arrival time and the amplitude of the wave and for the maximum wave run-up, where the friction terms play an important role. This test shows that the proposed hyperbolic strategy, the chosen breaking mechanism, and the standard friction term from classical shallow water equations perform adequately for the proposed hyperbolic system.

Solitary Wave Over a Gaussian Obstacle
This test run using the FG model is a solitary wave over a Gaussian obstacle in Ω = [−5, +35] × [−10, +10]. We consider the bathymetry given by where A g = 0.1 and σ g = 1. The initial soliton is taken from Test 4. has been carried out by considering a refined mesh, M2, made of 400 × 200 elements, using again N = 3. The 1D cut for the water depth at t = 12 is reported in Fig. 13. The numerical solutions obtained on both meshes match perfectly, hence mesh convergence has been reached for this test. Figure 14 shows the snapshots of the free-surface A = h + b at different times, t ∈ {0, 2, 5, 12}. As expected, the soliton is propagated growing in amplitude in correspondence to the obstacle and generating a set of transmissive waves behind the main wave.

Solitary Wave Impinging on a Conical Island
A set of 2D experiments with solitary waves was carried out at the Coastal and Hydraulics Laboratory, Engineer Research and Development Center of the U.S. Army Corps of Engineers, [16]. The laboratory experiment consisted on an idealized representation of Babi Island in the Flores Sea, Indonesia. It should be mentioned that the waves generated in the laboratory are dispersive; hence, as many authors have already done, [51,74,105], it constitutes an almost ideal problem to test the accuracy of the model.
As done in Sect. 4.1.6, the initial condition for ρ, v i , w and η, is a solitary wave of amplitude A = 0.06, H = 0.32 computed for the original SGN model and centred at x sol = 0. The wave propagates until t = 15 and the breaking mechanism is automatically activated with the same parameters than in all the previous numerical test, b 1 = 0.5, b 2 = 0.25. The relaxation speed λ is set to 300. Since the model is constructed with smooth concrete, the friction with the topography is not considered here. Finally, free-outflow boundary conditions are considered and the CFL number is set to 0.9.
Similarly to the test in Sect. 4.1.6, the a posteriori subcell finite volume limiter is employed in the presence of spurious oscillations as well as for values of the total water depth smaller than 10 −1 . During the entire simulation, the subcell finite volume limiter was active for less than the 4 % of total elements of the computational domain. That ensures using a robust second-order scheme in a tiny part of the computational domain.
As expected, the approximated solution, Fig. 15, shows two wavefronts splitting in front of the island and colliding behind it. In Fig. 16, it can be found that, overall, the simulated water surface fluctuations agree well with the measured data for both the maximum amplitudes and the arrival time of the waves. Therefore, the proposed methodology seems to be able to reproduce measured wave propagation over an uneven 3D bottom accurately.

Numerical Convergence Study and Evolution of Curl Errors
In this section we use the stationary radial solution found in Sect. 2.2.2 in order to carry out a numerical convergence study of the high order ADER schemes proposed in this paper and in order to measure the errors in the involution constraint ∇ × p as a function of the GLM curl cleaning speed a c .
The  Table 4, together with the observed order of accuracy. From the obtained results we can conclude that, on sufficiently fine meshes, our high order ADER-DG schemes achieve a numerical convergence order of N to N + 1 for all quantities.
We now run this test problem again for a longer time until t = 2 using an ADER-DG scheme with N = 4 on a uniform Cartesian mesh of N x = N y = 32 elements and report the time series of the L 2 error in the curl constraint ∇ × p for the original system (15)- (19) and of the augmented GLM system (35)-(40) as a function of the GLM cleaning speed a c . From the results shown in Fig. 17 we can deduce that the curl errors decrease with increasing cleaning speed, as expected. A set of radial cuts through this stationary vortex-type problem is shown in Fig. 18, where the numerical solution obtained with the high order ADER-DG scheme is compared against the exact solution of the problem at time t = 2. We can note an excellent agreement between numerical and exact solution.

Gray Soliton
In this section we solve the gray soliton test problem, which was introduced and solved with second order TVD finite volume schemes in [38]. In this paper we employ very high order ADER-DG schemes instead. The initial condition is given by and  According to [38] the parameters of the test problem are chosen as follows: b 1 = 1.5, b 3 = 1, U = 2, while the parameters of the hyperbolic model are set to β = 10 −3 , λ = 2000 and a c = 0, since the test problem is one-dimensional. Simulations are run until a final time of t = 20, when the solitary wave has traveled one period through the domain and has returned to its initial location, so that the exact solution is given by the initial condition. A comparison of the high order ADER-DG solution and the exact solution is shown in Fig. 19. For all depicted variables, which are representative for the system, we can note an excellent agreement between the numerical solution and the exact one. Note that the nominally fifth order ADER-DG scheme with N = 4 used in this paper needed only 5,400 degrees of freedom to resolve the solitary wave in x direction, while for the second order TVD scheme employed in [38] a very fine mesh composed of 200,000 elements in x direction was needed.
The following L 2 errors were measured at the final time t = 20: for the density, ρ, the error was 6.0462 · 10 −3 , for the first velocity component, v 1 , it was 6.3848 · 10 −3 , for η it was 6.0337 · 10 −3 and for p 1 it was 7.3374 · 10 −3 . These results clearly illustrate the benefit of very high order DG schemes for the numerical solution of hyperbolic reformulations of nonlinear dispersive systems over classical second order TVD finite volume schemes.

Dispersive Shock in One Space Dimension
In order to verify our proposed numerical approach and its implementation also in the presence of dispersive shock waves, we first compare the numerical solution obtained with a fourth order ADER-DG scheme (N = 3) in the domain Ω = [−250, +250]×[0, +25] on a uniform Cartesian mesh composed of 10080 × 4 elements with the exact solution of the Whitham modulation equations that describe the temporal evolution of the envelope of the dispersive shock wave, see [38]. The parameters of the hyperbolic model are set to β = 2·10 −5 , λ = 300 and a c = 0, since the test problem is one-dimensional. In Fig. 20 we show a comparison of numerical and exact solution for the fluid density ρ and for the velocity component v 1 at a final time of t = 70. For both variables, an excellent agreement between the exact and the numerical solution can be observed for the envelope.

Dispersive Shock in Two Space Dimensions
Last but not least, we also carry out a qualitative comparison against experimental results in two dimensions. We start from a fluid initially at rest with v k = 0, p k = ∂ k ρ, w = 0 and with ρ L = 2, ρ R = 1, R = 3 and the smoothing parameter σ = 0.05. The computational domain Ω = [−8, +8] 2 is discretized with 256 × 256 uniform Cartesian elements with polynomial approximation degree N = 3 and simulations are run until a final time of t = 0.24. Also for this test case the parameters of the hyperbolic model with GLM curl cleaning are set to β = 2 · 10 −5 , λ = 300 and a c = 20. In Fig. 21 we compare our numerical results against experimental observations published in [65] for a two-dimensional dispersive shock wave developing in a quantum fluid (Bose-Einstein Condensate-BEC). We can observe a good qualitative match between the two solutions.

Conclusion
In this paper we have focused on the solution of first order hyperbolic reformulations of nonlinear dispersive systems using a high order fully discrete one-step ADER discontiuous Galerkin methodology. Two different models derived from an extended Lagrangian variational principle have been considered. The hyperbolic reformulation of the Serre-Green-Naghdi model with flat bottom [56] has been extended to the 2D framework with variable bottom topography. The linear dispersion analysis reported in the present paper allows an easy comparison with other hyperbolic models of nonlinear dispersive water waves. As second system we have selected the defocusing nonlinear Schrödinger equation rewritten in the hyperbolic reformulation proposed in [38]. The extension of the original 1D model to the multi-dimensional framework has shown the necessity of including a new term to restore the Galilean invariance of the system. Moreover, in this paper the model has also been augmented by a hyperbolic GLM curl cleaning approach so that we arrive at a structure preserving scheme that reduces the curl errors that are produced by the underlying numerical method. Let us note that, as an additional benefit of including the curl cleaning terms, the augmented system becomes strongly hyperbolic when the second and third components of the vector field p are non zero, reverting to the eigenstructure of the original system in [38], which was only weakly hyperbolic in the multi-dimensional case in the absence of the GLM curl cleaning terms. The path-conservative ADER-DG methodology has been briefly recalled and combined with an a posteriori subcell finite volume limiter to avoid spurious oscillations in the numerical solution in the presence of discontinuities or strong gradients. For the hyperbolic model of nonlinear dispersive water waves (1)-(5) the well-balanced property of the employed path-conservative scheme has been rigorously proven for the case N = 0 for arbitrary bottom topographies. A careful validation of the schemes has been carried out for both systems of equations showing numerical convergence tables of up to sixth order of accuracy.
A set of numerical test has been run for each system including a comparison with exact and numerical reference solutions, as well as with available experimental data. The obtained results confirm the effectiveness of using very high order accurate numerical schemes against classical low order finite volume methods which, even in the one-dimensional case, require a very fine grid leading to a prohibitive computational cost in higher space dimensions.
In the future, we plan to apply the numerical method developed in this paper also to other nonlinear dispersive systems, such as the Navier-Stokes-Korteweg system, for which a hyperbolic reformulation is possible by combining the GPR model of continuum mechanics [17,47,88] with the hyperbolic reformulation of the nonlinear Schrödinger equation introduced in [38] and studied in the present paper. As further research, we also plan to develop novel structure-preserving semi-implicit schemes on staggered meshes to solve the systems under investigation in this paper, following the ideas outlined in [9,11,18,19], with the aim to preserve the curl constraint exactly at the discrete level. Funding Open access funding provided by Università degli Studi di Trento within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. not influence the corresponding Euler-Lagrange equations. The incompressibility equation and kinematic boundary conditions imply the mass conservation law in the form h t + div(hv) = 0, hv = b+h b vdz. (97) In the following, 'dot' means the material derivative along the average velocity. For example, Under additional assumption -the vorticity is of order ε β , with β > 1-one can easily obtain (see [2] for details) Then, up to ε 2 order terms one has [2]: It allows us to write the Lagrangian in the form (we return back to the dimension variables): The corresponding Hamilton's action between the time instant t 0 and t 1 is then The Euler-Lagrange equations for (101) are obtained under the mass conservation constraint (97). Following [56], the extended Lagrangian is defined aŝ Here λ is a large constant having the dimension of squared velocity. The corresponding variational technique to find the Euler-Lagrange equations can be found in [36], [38]. The obtained governing equations are yet quite complicated for the numerical study and contain non-conservative products. We simplify now the Lagrangian (102) under the hypothesis that Let us now take Then Hence, it is sufficient to take A = 1 2 a c ψ × p to have the energy conservation. Finally, in the thermodynamically compatible GLM cleaning procedure we must take , a c > 0.
In the main text we present this model also in coordinate form.