High order ADER-DG schemes for the simulation of linear seismic waves induced by nonlinear dispersive free-surface water waves
HHigh order ADER-DG schemes for the simulation of linearseismic waves induced by nonlinear dispersive free-surface waterwaves
C. Bassi ( b ) , S. Busto ( a ) , M. Dumbser ( a )1 ( a ) Laboratory of Applied Mathematics, DICAM, University of Trento, via Mesiano 77, 38123 Trento, Italy ( b ) MOX–Modelling and Scientific Computing, Dipartimento di Matematica, Politecnico di Milano, PiazzaLeonardo da Vinci 32, 20133 Milano, Italy
Abstract
In this paper, we propose a unified and high order accurate fully-discrete one-step ADER Discontin-uous Galerkin method for the simulation of linear seismic waves in the sea bottom that are generatedby the propagation of free surface water waves. In particular, a hyperbolic reformulation of the Serre-Green-Naghdi model for nonlinear dispersive free surface flows is coupled with a first order velocity-stressformulation for linear elastic wave propagation in the sea bottom. To this end, Cartesian non-conformingmeshes are defined in the solid and fluid domains and the coupling is achieved by an appropriate time-dependent pressure boundary condition in the three-dimensional domain for the elastic wave propagation,where the pressure is a combination of hydrostatic and non-hydrostatic pressure in the water column abovethe sea bottom. The use of a first order hyperbolic reformulation of the nonlinear dispersive free surfaceflow model leads to a straightforward coupling with the linear seismic wave equations, which are alsowritten in first order hyperbolic form. It furthermore allows the use of explicit time integrators with arather generous CFL-type time step restriction associated with the dispersive water waves, compared tonumerical schemes applied to classical dispersive models that contain higher order derivatives and typi-cally require implicit solvers. Since the two systems that describe the seismic waves and the free surfacewater waves are written in the same form of a first order hyperbolic system they can also be efficientlysolved in a unique numerical framework. In this paper we choose the family of arbitrary high orderaccurate discontinuous Galerkin finite element schemes, which have already shown to be suitable for thenumerical simulation of wave propagation problems. The developed methodology is carefully assessed byfirst considering several benchmarks for each system separately, i.e. in the framework of linear elasticityand non-hydrostatic free surface flows, showing a good agreement with exact and numerical referencesolutions. Finally, also coupled test cases are addressed. Throughout this paper we assume the elasticdeformations in the solid to be sufficiently small so that their influence on the free surface water wavescan be neglected.
Keywords: hyperbolic equations, ADER schemes, discontinuous Galerkin finite element methods, hy-perbolic reformulation of the Serre-Green-Naghdi model, linear elastic wave equations in velocity-stressformulation, coupling of nonlinear dispersive water waves with linear elastic waves
The physical phenomenon we are interested in is the generation of seismic waves in the sea bottom dueto the propagation of free surface water waves on the sea surface in near coastal regions. In view of thephysical characteristics of the two materials involved, the wave speeds in the solid medium and the freesurface wavespeed in the fluid can differ by up to two orders of magnitude. To simulate such a complexsituation, we propose the use of a high order accurate fully-discrete one-step ADER-DG scheme on nonconforming meshes, which solve a coupled set of two first order hyperbolic systems. The first model is ahyperbolic reformulation [54] of the Serre-Green-Naghdi model for the description of non-hydrostatic free- Corresponding author.Email addresses: [email protected] (C. BAssi), [email protected] (S. Busto), [email protected] (M.Dumbser) a r X i v : . [ m a t h . NA ] A p r urface flows, while the second model consists in a classical first order hyperbolic model of linear elasticityin velocity-stress formulation for the description of linear seismic waves propagation [111, 7, 75].The most straightforward way to model water wave propagation could be to consider a fully three-dimensional free surface flow model, like for example those developed in a series of papers by Casulli etal. [21, 25, 15, 22, 24, 23], where an efficient semi-implicit method for fully three-dimensional hydrostaticand non-hydrostatic free surface flows has been proposed. Since, however, we are interested in thepropagation of free surface sea waves near the coast, where the typical horizontal length scales are farlarger than the typical vertical length scales, also simplified shallow water-type equations (SWE) may be asuitable alternative, substantially reducing computational complexity and increasing the efficiency of thefinal methodology compared to fully three-dimensional models. A wide variety of phenomena associatedto the propagation of water waves can be successfully described employing the classical shallow waterequations [35]. Nevertheless, the SWE are unable to reproduce non-hydrostatic effects and propagation ofsolitary and dispersive waves. For this reason, we need to go beyond the classical shallow water equations,looking for more sophisticated dispersive systems. Starting from the pioneering work [11], in which afirst 1D Boussinesq-type model is derived under the assumptions of weak dispersion, weak non-linearityand flat bottom topography, different dispersive systems have been proposed in the literature. Amongthem we recall the Peregrine system, [88], where the model in [11] is extended to two dimensions and tonon flat bottom topography (still maintaining the weak nonlinearity hypothesis), the Serre model, [95],where a 1D fully non-linear approach is presented for flat bottom topography (still keeping the weakdispersion assumption), the Serre-Green-Naghdi (SGN) model, [65, 94], where an extension of [95] to twodimensions for arbitrary bottom is provided and [28], where a derivation of the model in [65] is built,using asymptotic expansion and irrotationality. Notice that, in addition to the models just mentioned,more advanced models with improved dispersion characteristics (see [79, 80, 85, 77, 78]) and models whichinclude additional physical phenomena with respect to the classical formulation for inviscid flows (see thereview [71]) have been also proposed in the literature.An important distinguishing feature of dispersive models is that they contain higher order space andspace-time derivatives. Consequently, their numerical discretization becomes particularly difficult anda severe time step restriction arises when explicit time integration schemes are employed. A possiblesolution to this major drawback is the introduction of augmented first order systems, such as in [14, 59].Also in this case, however, the hyperbolicity of the SW equations is lost, thus leading to the necessity ofsolving an elliptic equation for the dispersive part at each time-step. A completely different approach isinstead adopted in [54], where a hyperbolic approximation of the non-hydrostatic system in [14] and ofthe SGN model with mild bottom approximation are proposed. The hyperbolic approximation proves tohave a dispersion relation which is very similar to the dispersion relations associated to the original non-hyperbolic models and allows to obtain very accurate numerical results. At the same time hyperbolicityallows to easily implement the model in the context of high-order finite volume (FV) and discontinuousGalerkin (DG) schemes and to realize efficient numerical simulations also in multiple space dimensions.Due to these evident advantages, we have decided to employ this approach for the water waves part inour coupled simulations. Notice that the introduction of hyperbolic reformulations is not a novelty andcomes from the pioneering work by Cattaneo [26], where the second order terms in the heat equation arerewritten as relaxation terms, while a hyperbolic approximation has been proposed for the first time inthe context of dispersive systems in [58] (precisely for the Serre model [95]). Further developments of[54] are presented in [5], where an hyperbolic reformulation of SGN model for arbitrary bathymetry inintroduced. Moreover, in [53], an even more general formulation has been proposed, including furtherdispersive Boussinesq-type systems, as the models of Yamazaki et al., [112], and Peregrine, [88, 79],which are built neglecting some convective terms in the equation for the averaged vertical velocity andin an auxiliary equation accounting for the spatial variation of the mean horizontal velocity. For furtherfirst order hyperbolic reformulations of dispersive and dissipative systems, see also [37] and [89, 48].Concerning the propagation of seismic waves, we put ourselves in the context of linear elasticity andwe adopt a first order velocity-stress formulation, which has the advantage of being hyperbolic, see e.g.[75, 69].A key idea of the present paper is the mathematical description of both wave propagation problems,i.e. the linear seismic waves as well as the nonlinear dispersive water waves, at the aid of first order2yperbolic systems. For this reason, we can easily couple both models with each other and solve themnumerically in the time domain by employing high order accurate discontinuous Galerkin (DG) finiteelement methods in a straightforward manner. We decide, in particular, to use an explicit high orderADER-DG scheme [69], for which the usual CFL condition holds, i.e. with ∆ t proportional to ∆ x , thusavoiding higher powers of ∆ x that would be typical for explicit time discretizations of Boussinesq-typeequations with higher order spatial and temporal derivatives.The DG method has been introduced for the first time in [90], for the solution of a neutron transportequation. It has been then extended to time-dependent, multidimensional and nonlinear hyperbolic prob-lems in [27, 33, 30, 29, 31]. Besides, in [6, 32] the Local Discontinuous Galerkin (LDG method) has beenproposed in order to solve convection-diffusion equations. This approach involves rewriting the secondorder equations as an augmented non-hyperbolic first order system and the subsequent discretization ofthis augmented system with the DG method. The DG method has been applied for the first time toequations containing higher order derivatives in [114], where the LDG method is used for the resolutionof linear dispersive Korteveg-de-Vries (KdV) equations, containing up to third order spatial derivatives.Extensions to linear equations with derivatives up to fifth order and nonlinear dispersive equations are in-stead presented in [113, 76]. Applications of the DG method to the solution of nonlinear Boussinesq-typedispersive equations have been introduced in [56, 55, 52]. Notice that the severe time step restrictionsdue to higher-order spatial derivatives are overcome in [42], where a fully implicit space-time DG methodis applied to both linear third order KdV equations and nonlinear Boussinesq-type systems.High order of accuracy in space is straightforward to obtain in the DG framework, while attaininghigh order in time is still a very active field of research. A successful approach consists in using thealready mentioned space-time DG methods, [92, 93, 109, 110, 97, 98, 72, 99, 18]. An alternative, thatwill be followed in this work, are the ADER-DG schemes, first put forward in [41] in the context of FVmethods and generalized in [39] to the unified P N P M framework for arbitrary high order accurate one-step FV and DG schemes. Classical ADER methods (Arbitrary high order DErivative Riemann problem)have been proposed by Millington et al., [82], and Toro et al., [104], in the framework of finite volumemethods. The methodology is based on the resolution, at the cell interfaces, of a generalized Riemannproblem with piecewise polynomial initial conditions, built using a nonlinear reconstruction (e.g. ENO orWENO methods) that circumvents Godunov’s theorem. Then, space-time integration on an appropriatecontrol volume is performed using a Taylor series expansion in time, where time derivatives are replacedby spatial derivatives following the Cauchy-Kovalevskaya procedure. Further developments of classicalADER schemes, including their extension to the DG framework, can be found in [106, 102, 47, 46, 63,19, 17, 36, 9] and references therein. The main inconvenient of classical ADER methods is the needof the cumbersome Cauchy-Kovalevskaya procedure. The novel ADER-DG methodology presented in[41, 39] avoids that step by employing a new element-local space-time DG predictor, which leads tomore efficient algorithms. Since then, ADER-DG has been used to solve many different models, such ascompressible flows [38], hyperbolic systems in general relativity [44, 57, 43] or non conservative hyperbolicsystems for geophysical flows [40]. Moreover, the ADER-DG method has also proven to be well suitedfor the simulation of seismic waves problems, see [69, 45, 70, 34, 46]. For other high order discontinuousGalerkin finite element schemes applied to elastic wave equations, see e.g. [66, 2, 1, 3, 4]. Finally, in [54]the ADER-DG method has also been successfully applied to the solution of hyperbolic reformulations ofdispersive models. Due to these considerations, the ADER-DG method appears to be a suitable choicefor the discretization of the coupled system proposed in the present work.The paper is organized as follows. In Section 2, the mathematical models employed for the descriptionof both seismic and free surface water waves are recalled. Moreover, the boundary conditions to be seton the interface between the fluid and solid domains are defined. Section 3 is devoted to the descriptionof the high order one-step ADER-DG scheme on Cartesian grids. In Section 4, the numerical methodis validated. First, some classical benchmarks are solved independently for each of the two systems ofequations. Then, two test cases of the coupled problem are presented. Finally, in Section 5, we draftsome conclusions and perspectives. 3 Governing equations
As already mentioned within the introduction, we aim at simulating the effect of surface water waveson the generation and propagation of seismic waves on the sea bottom. To this end, we couple the twosystems that are recalled in this section: a non-hydrostatic dispersive shallow water-type model for thepropagation of the free surface water waves and a linear elasticity model for the seismic waves propagatingin the solid domain below the sea floor.
During the last decades many non-hydrostatic models aiming at characterising non-hydrostatic free sur-face flows have been successfully developed. Taking into account the features of the flow to be modelled,we will employ the Serre-Green-Naghdi (SGN), [64, 65] model, and a dispersive model recently proposedby Sainte-Marie et al. (SM) in [14]. Both models are used in combination with the so-called mild bot-tom approximation. It is important to remark that both systems can be rewritten in a one-parameterdependent unified formulation, [54], as ∂h∂t + ∇ · ( h u ) = 0 , (1) ∂h u ∂t + ∇ · ( h u ⊗ u ) + ∇ · (cid:18) gh I + hp I (cid:19) + ( gh + γp ) ∇ z b = 0 , (2) ∂hw∂t + ∇ · ( h u w ) = γp, (3) w + 12 h ∇ · u − u · ∇ z b = 0 , (4)where we have denoted h the water depth, u = ( u x , u y ) the horizontal velocity vector of the fluid, w theauxiliary variable of the averaged vertical flow velocity, p the depth-averaged non-hydrostatic correctionfor the pressure, g = 9 .
81 the gravity acceleration and z b the vertical coordinate of the bottom bathymetry.Moreover, the gradient and divergence operators considered refer to the horizontal plane, ∇ = ( ∂ x , ∂ y ) T ,neglecting the vertical variable and I is the identity matrix of dimension two. Substituting γ = , theformer system would provide the SGN model whereas setting γ = 2 leads to the SM model.Following [54], the first order unified reformulation of the SGN and SM models reads ∂h∂t + ∇ · ( h u ) = 0 , (5) ∂h u ∂t + ∇ · ( h u ⊗ u ) + ∇ · (cid:18) gh I + hp I (cid:19) + ( gh + γp ) ∇ z b = 0 , (6) ∂hw∂t + ∇ · ( h u w ) = γp, (7) ∂hp∂t + ∇ · ( h u p ) + c (2 w + h ∇ · u − u · ∇ z b ) = 0 . (8)with c the artificial sound speed, c = α √ gH , H the average still water depth and α >
0. The system(5)-(8) represents the conservation of mass and momentum and is furthermore augmented by two PDEsfor the auxiliary variables p and w , with appropriate source terms, that allow the system to relax, for c → ∞ , towards the original SGN or SM models, depending on the value of the parameter γ . Moreover,assuming steady bathymetry, the system can be completed with the following equation: ∂ t z b = 0 . (9)It is important to note that the above system is depth averaged and, thus, the vector of spatial coordinates x ∈ Ω w ⊂ R , where Ω w is the two-dimensional computational domain used for the simulation of thewater wave propagation (see Figure 1). The fact that the equations are depth averaged, jointly with theuse of a two-dimensional domain, makes this model particularly interesting concerning its computationalefficiency, compared to a more complete three-dimensional non-hydrostatic formulation.4igure 1: Two-dimensional computational domain for the simulation of the water wavepropagation. h denotes the water depth, H is the still water depth, A = h − H corre-sponds to the surface elevation with respect to the still water depth, η is the total surfaceelevation, and z b gives the bottom bathymetry. Assuming the sea bottom to be a homogeneous isotropic elastic material under small deformations, itcan be modelled using the linear elasticity equations [67, 8, 7, 75]. The classical formulation of linearelasticity is a second order vector wave equation, but it can be rewritten as a first order hyperbolic systemin velocity-stress formulation, which leads to an easier and more direct coupling with the model governingthe dispersive water waves shown in the previous section. Accordingly, the final hyperbolic system forthe linear elastic wave propagation reads ∂ σ ∂t − λ ( ∇ · v ) I − µ (cid:0) ∇ v + ∇ v T (cid:1) = 0 , (10) ∂ρ v ∂t − ∇ · σ = 0 , (11)where the first equation, (10), is the Hooke law expressed in terms of the two Lam´e constants λ and µ andthe second equation, (11), represents the conservation of momentum. Here, the time is again denoted by t and the spatial coordinate is x ∈ Ω e ⊂ R , with Ω e a 3D computational domain used for the simulationof the seismic wave propagation. Furthermore, in the above system, σ = σ T is the symmetric stresstensor, ˜ σ = ∂ σ ∂t corresponds to the linearized part of the first Piola-Kirchhoff stress tensor, u is the vectorof the velocity field and ρ is the density. To simulate the effect of free surface waves on the generation and propagation of seismic waves on thesea bottom, the non-hydrostatic and the linear elastic wave models need to be coupled. Throughout thispaper we assume that the coupling is only done in a one-way manner, i.e. the free surface waves lead to atime-dependent pressure on the sea floor, which generates low frequency and low amplitude linear elasticwaves in the solid medium below the sea bottom. Instead, we assume that the elastic deformations of thesolid are very small and therefore do not couple back to the surface water waves via a time-dependentbottom geometry. The non-hydrostatic description of the free surface water waves propagation thereforeprovides appropriate boundary conditions for the stress tensor on the upper boundary of the domain Ω e ,which will be denoted by ∂ Ω e,u in the following. In particular, we assume zero shear stress and continuityof the normal stress on ∂ Ω e,u , σ xz = σ yz = 0 , σ zz = ρgh. (12)As already stated before, due to their very low amplitude, the linear elastic waves are not coupledback to the non-hydrostatic shallow water solver, although one might imagine a fully coupled systemby considering, in (5)-(9), ∂ t z b = v · n z , where n z is the unit normal vector in z direction and v is thevelocity of the solid. 5 .4 Unified writing of the models To provide a unified description of the numerical scheme for the two models considered, we rewrite thesystems of equations to be solved in the general form ∂ t Q + ∇ · F ( Q ) + B ( Q ) · ∇ Q = S ( Q ) , (13)where Q = Q ( x , t ) is the vector of unknowns, F ( Q ) is the nonlinear flux tensor, B ( Q ) · ∇ Q is agenuinely non-conservative term, and S ( Q ) is an algebraic source term. Therefore, for the hyperbolicnon-hydrostatic model (5)-(9), we have Q = hhu hu hwhpz b , F ( Q ) = hu hu hu + gh + hp hu u hu u hu + gh + hphwu hwu hu (cid:0) p + c (cid:1) hu (cid:0) p + c (cid:1) , B ( Q ) · ∇ Q = gh + γp ) ∂ x z b ( gh + γp ) ∂ y z b c [ − u ∂ x h − u ∂ y h − u ∂ x z b + u ∂ y z b )]0 , S ( Q ) = γp − c w . (14)On the other hand, the definition of Q = ( σ xx , σ yy , σ zz , σ xy , σ yz , σ xz , u , u , u ) T and B ( Q ) · ∇ Q = − ( λ + 2 µ ) ∂ x u − λ∂ y u − λ∂ z u − λ∂ x u − ( λ + 2 µ ) ∂ y u − λ∂ z u − λ∂ x u − λ∂ y u − ( λ + 2 µ ) ∂ z u − µ∂ x u − µ∂ y u − µ∂ y u − µ∂ z u − µ∂ x u − µ∂ z u − ρ ( ∂ x σ xx + ∂ y σ xy + ∂ z σ xz ) − ρ ( ∂ x σ xy + ∂ y σ yy + ∂ z σ yz ) − ρ ( ∂ x σ xz + ∂ y σ yz + ∂ z σ zz ) (15)leads to the linear elasticity model (10)-(11). 6 Numerical discretization
The high order accurate fully-discrete one-step ADER discontinuous Galerkin methodology (ADER-DG),[39, 50], is used in order to discretize the two models considered, namely the hyperbolic reformulation ofthe SGN equations (HSGN), presented in Section 2.1, and the linear elasticity system recalled in Section2.2. ADER-DG methods fall into the framework of the general P N P M schemes proposed in [39], thatextend the local predictor ADER methodology presented in [41] for FV also to the DG framework. Moreprecisely, we focus on the pure DG case, where N = M , which has shown to be appropriate to solvelinear and non-linear hyperbolic conservation laws. In this section, we provide a brief summary of themethod on Cartesian grids. Further developments of this methodology, including the use of unstructuredmesh and Cartesian grids with adaptative mesh refinement (AMR) employed to solve a wide variety ofhyperbolic PDEs, in both the Eulerian and the Lagrangian framework, can be found, for instance, in[38, 49, 50, 115, 51, 60, 44, 16] and references therein.Before describing the numerical method to be employed, we define a discretization of the compu-tational domain Ω using a Cartesian grid made of elements of the form Ω i = (cid:2) x i − ∆ x, x i + ∆ x (cid:3) × (cid:2) y i − ∆ y, y i + ∆ y (cid:3) × (cid:2) z i − ∆ z, z i + ∆ z (cid:3) , with x i = ( x i , y i , z i ) the barycentre of cell Ω i and ∆ x ,∆ y , ∆ z the cell size on each spatial coordinate direction. Next, following the classical DG approach, weassume that the space of discrete solutions of (13) is generated by spatial basis functions φ k constructedas the tensor product of piecewise polynomials up to degree N . In particular, we consider the orthogonalLagrange interpolation polynomials passing through the Gauss-Legendre quadrature points of a N + 1Gauss quadrature formula. Then, within each element Ω i the discrete solution of the system can bewritten as u h ( x , t n ) = φ l ( x ) ˆ u nl , x ∈ Ω i , (16)where the classical Einstein summation convection is employed, ˆ u nl denote the degrees of freedom of thesolution, and l is a multidimensional index referring to the one-dimensional basis functions, φ l m , on areference element Ω ref = [0 , ≤ ξ, η, ζ ≤ x = x i − ∆ x + ξ ∆ x , y = y i − ∆ y + η ∆ y , and z = z i − ∆ z + ζ ∆ z ,respectively. Since the chosen basis functions are not time dependent, the direct use of a classical DGapproach would result in a low order scheme in time. To attain high order of accuracy in time, we usethe ADER-DG methodology which can be divided into two main steps: • Local space-time predictor computation. System (13) is solved “in the small” using a locally im-plicit space-time discontinuous Galerkin scheme on each element, which neglects iterations betweenneighbour cells. • Explicit update of the solution using a one-step corrector. The space-time predictor is used intothe weak formulation of (13) which takes into account the fluxes between cells and provides thesolution of the system at the new time instant.We come now to further detail each step.
To determine the local space-time predictor solution, q h ( x , t ), which will lead to a high order scheme inspace and time avoiding the cumbersome Cauchy-Kovalewskaya procedure used in the original ADERmethods [105, 107, 101, 108], we employ the weak formulation in space-time proposed in [41, 39]. Let usconsider the space-time test functions, θ k = θ k ( x , t ), built as the product of the already introduced nodalspatial basis functions and an additional one dimensional basis function for the time dependency. withthe additional transformation for the reference time τ given by t = t n + τ ∆ t . Then, multiplying (13) by θ k and integrating over the space-time control volume, Ω i × (cid:2) t n , t n +1 (cid:3) , yields t n +1 (cid:90) t n (cid:90) Ω i θ k ∂ t q h d x dt + t n +1 (cid:90) t n (cid:90) Ω i θ k ∇ · F ( q h ) d x dt + t n +1 (cid:90) t n (cid:90) Ω ◦ i θ k B ( q h ) · ∇ q h d x dt = (cid:90) Ω i θ k S ( q h ) d x dt. (17)7ithin an implicit space-time DG method, [109, 110], the weak formulation (17) would be now integratedby parts in space and time to provide the solution at the new step. However, we are just interested inobtaining a local approximation of the predictor so we only integrate by parts in time, neglecting theinteraction between neighbours, (cid:90) Ω i θ k ( x , t n +1 ) q h ( x , t n +1 ) d x − (cid:90) Ω i θ k ( x , t n ) u h ( x , t n ) d x − t n +1 (cid:90) t n (cid:90) Ω i ∂ t θ k q h d x dt + t n +1 (cid:90) t n (cid:90) Ω i θ k ∇ · F ( q h ) d x dt + t n +1 (cid:90) t n (cid:90) Ω ◦ i θ k B ( q h ) · ∇ q h d x dt = (cid:90) Ω i θ k S ( q h ) d x dt. (18)Moreover, in (18), we have taken into account that the predictor at time t n is given by the degrees offreedom of the solution at the previous time step, u h ( x , t n ), thus respecting the causality principle. Thisfact can be seen as the realization of an upwinding approximation in time. Therefore, the above nonlinearsystem has as only unknown the degrees of freedom of the space-time expansion, q h ( x , t ) = θ k ( x , t ) ˆ q k (19)and can be solved locally at each cell using a discrete Picard iteration procedure. Since the Picarditeration matrix is nilpotent, it will converge in at most N + 1 iterations, as it has been proven in [68] forhomogeneous linear conservation laws. The convergence of the Picard iteration for nonlinear systems ofconservation laws was proven in [16]. The solution of (18) constitutes the only (element-local) implicitstep on the whole ADER-DG algorithm. The solution q h , obtained at the predictor step, does not account for the neighbouring flux contributions,so it can not be used as the solution of the PDE system at time t n +1 . To correct this issue we employ anexplicit one-step DG approach. We first multiply the governing PDE system (13), by the test functions φ k and we then integrate over a space-time control volume Ω i × [ t n , t n +1 ], obtaining the following weakproblem t n +1 (cid:90) t n (cid:90) Ω i φ k ( ∂ t Q + ∇ · F ( Q ) + B ( Q ) · ∇ Q ) d x dt = t n +1 (cid:90) t n (cid:90) Ω i φ k S ( Q ) . (20)Taking into account (16), integrating the flux divergence term by parts in space and the time derivativeby parts in time yields (cid:90) Ω i φ k φ l d x (cid:0) ˆ u n +1 l − ˆ u nl (cid:1) + t n +1 (cid:90) t n (cid:90) ∂ Ω i φ k G (cid:0) q − h , q + h (cid:1) · n dS dt + t n +1 (cid:90) t n (cid:90) ∂ Ω i φ k D (cid:0) q − h , q + h (cid:1) · n dS dt − t n +1 (cid:90) t n (cid:90) Ω i ∇ φ k · F ( q h ) d x dt + t n +1 (cid:90) t n (cid:90) Ω ◦ i φ k B ( q h ) · ∇ q h d x dt = (cid:90) Ω i φ k S ( q h ) d x dt, (21)where n denotes the outward unit normal at the cell boundary, ∂ Ω i , and q h = q h ( x , t ) is the local space-time predictor already introduced in Section 3.1. Since we are using a discontinuous Galerkin scheme, thebasis functions are allowed to jump across cell interfaces. To account for the discontinuities arising in thesecond term of (21), we make use of a Riemann solver at the element interfaces, see e.g. [103]. The initialcondition for the numerical flux function is then given by the right and left states q − h , q + h computed at8he predictor step which yields the order of accuracy sought. In particular, within this work, we considerthe Rusanov numerical flux function, hence G (cid:0) q − h , q + h (cid:1) · n = 12 (cid:0) F ( q + h ) + F ( q − h ) (cid:1) · n − s max (cid:0) q + h − q − h (cid:1) . (22)Furthermore, we also need to develop a proper discretization of the non conservative products at theelement boundaries. To this end, we consider the works [87, 20, 86, 84], based on the theory of Dal Maso,Le Floch and Murat [81] on nonconservative hyperbolic PDE systems, and their extension to higher orderDG schemes in [91, 40]. Within this framework it is usual to build also so-called well-balanced schemes[10, 62, 74]. Accordingly, the third term in (21) is approximated with a path integral in phase-spacebetween the two extrapolated values related to the face, q − h and q + h , D (cid:0) q − h , q + h (cid:1) · n = 12 (cid:90) B (cid:0) ψ ( q − h , q + h , s ) (cid:1) · n ds · (cid:0) q + h − q − h (cid:1) , (23)where we have used the linear segment path ψ = ψ ( q − h , q + h , s ) = q − h + s (cid:0) q + h − q − h (cid:1) , s ∈ [0 , . (24)As it can be seen in [83, 61], different paths could have been chosen to perform the former integralattending to special features of non conservative and source terms. Nevertheless, for the systems addressedin this work, the easiest straight line path already provides good results.From the computational point of view, it is important to remark that, as a consequence of thenodal tensor-product basis employed, the scheme can be written in a dimension by dimension fashionand integral operators are decomposed as the product of one-dimensional operators [16]. The resultingexplicit one-step ADER-DG scheme is conditionally stable with stability condition∆ t < CFL h min (2 N + 1) | λ max | , (25)with CFL < /d , where | λ max | is the maximum eigenvalue of the system, h min denotes the minimumcharacteristic mesh size and d is the number of space dimensions. The final time step for the coupledproblem is the lowest one between the timesteps associated to the two considered models. A major challenge concerning the coupling of the two models is the large discrepancy, in the spatiallength scales, between the elastic waves, with a typical propagation speed of about 3000 − m in water. To address this problem, we employ two non-conforming meshes withdifferent spatial resolutions. An initial 3D mesh made of rectangular bricks is first designed to cover Ω e .Then, on Ω w , a much finer grid is built as a refinement of the 2D mesh made up by the faces of the meshdesigned in Ω e , lying on the boundary ∂ Ω e,w (see Figure 2). Interpolation of the discrete solution fromΩ w onto the boundary ∂ Ω e,w is carried out using appropriate high order Gaussian quadrature formulas,[96], so that the boundary condition (12) can be imposed on the solid domain.Further boundary conditions for the solid media include periodic and free surface boundaries. Re-garding the last ones, the exact Riemann solver of Godunov can be employed, [69, 45, 100]. Accordingly,the opposite value for the incoming normal stress to the boundary is set. Finally, on the fluid domain weconsider either periodic or Dirichlet boundary conditions. This section is devoted to the assessment of the developed methodology. Initially, we address the twomathematical models independently, analysing the solution obtained for classical benchmarks of linear9igure 2: Sketch of the nonconforming mesh at the interface, ∂ Ω e,w . The 2D cells (grey)filling the fluid domain, Ω w , are built as a subgrid of the faces of the 3D elements dis-cretizing the solid domain, Ω e .elasticity and non-hydrostatic flows. Once the numerical method is validated, we present two showcasesof the coupled problem, reporting the seismic wave propagation generated by a soliton and by a sinusoidalwave train on the water surface. Linear elasticity is a well established research field, so many numerical tests can be found to validatethe proposed methodology. In what follows, we will first validate the numerical method using a p-s-wavetest, whose exact periodic solution is known. Then, we focus on a classical benchmark of seismic wavepropagation, the so-called Lamb problem, and on a stiff inclusion test, which accounts for large materialparameter variations. A final wave propagation test in three dimensions is also included.
Following [69, 99], a p − and s − wave test case is employed to verify the accuracy of the ADER-DGscheme. We consider the computational domain Ω = [ − . , . × [ − . , . ∈ R with periodic boundaryconditions in x and y directions and we define the initial condition Q ( x ,
0) = α r p sin (2 π n · x ) + α r s sin (2 π n · x ) , (26)with α = 0 . n = ( n x , n y ) = (1 , r p , r s the eigenvectors associated to the p − and s − waves, r p = (cid:0) λ + 2 µn x , λ + 2 µn y , µn x n y , − c p n x , − c p n y (cid:1) , r s = (cid:0) − µn x n y , µn x n y , µ (cid:0) n x − n y (cid:1) , c s n y , − c s n x (cid:1) , (27) c p = (cid:112) ( λ + 2 µ ) /ρ the p − wave speed and c s = (cid:112) µ/ρ the s − wave speed. Setting the material parametersto λ = 2, µ = 1, ρ = 1, leads to propagation velocities, c p = 2 and c s = 1. The former initial conditiongenerates a sinusoidal p − wave travelling along the diagonal direction, n = (1 , s − wave moving in the opposite direction. Taking t end = 3 √
2, the solution coincides with the initial datumand a convergence analysis can be performed. Table 1 shows the L errors, (cid:15) L , and the convergenceorder, O , obtained for N ∈ { , , , , } , M = N . The spatial grids were built using the same number ofelements in x and y directions, N y = N x , and the time step is computed at each iteration according tothe CFL condition. All variables attain the optimal order of convergence sought. The Lamb’s problem is a well known benchmark used to test numerical methods for linear elastic wavesand has first been put forward in [73]. In this paper we consider one of its variants, already analysed in[69, 48], where a rectangular domain Ω = [ − , × [ − ,
0] and a point source of the form f v ( x , t ) = ρ s a (cid:18)
12 + a ( t − t D ) (cid:19) exp (cid:16) a ( t − t D ) (cid:17) δ ( x − x s ) e y (28)10 N x (cid:15) L ( u ) (cid:15) L ( v ) (cid:15) L ( σ xx ) (cid:15) L ( σ yy ) (cid:15) L ( σ xy ) O ( u ) O ( v ) O ( σ xx ) O ( σ yy ) O ( σ xy ) Teor.3 10 1 . E −
02 1 . E −
02 2 . E −
02 3 . E −
02 1 . E −
02 - - - - - 415 2 . E −
03 3 . E −
03 5 . E −
03 6 . E −
03 1 . E −
03 4 .
10 3 .
99 4 .
08 4 .
16 4 . . E −
04 1 . E −
03 1 . E −
03 2 . E −
03 6 . E −
04 4 .
03 3 .
92 4 .
00 4 .
04 4 . . E −
04 2 . E −
04 3 . E −
04 4 . E −
04 1 . E −
04 4 .
02 3 .
95 3 .
98 4 .
00 4 . . E −
05 6 . E −
05 1 . E −
04 1 . E −
04 3 . E −
05 4 .
01 3 .
97 3 .
94 3 .
97 4 .
014 10 7 . E −
04 1 . E −
03 1 . E −
03 2 . E −
03 6 . E −
04 - - - - - 515 1 . E −
04 1 . E −
04 2 . E −
04 3 . E −
04 9 . E −
05 4 .
81 4 .
99 4 .
79 4 .
80 4 . . E −
05 3 . E −
05 6 . E −
05 7 . E −
05 2 . E −
05 4 .
90 5 .
01 4 .
85 4 .
87 4 . . E −
06 4 . E −
06 8 . E −
06 1 . E −
05 3 . E −
06 4 .
94 5 .
01 4 .
87 4 .
89 4 . . E −
07 1 . E −
06 2 . E −
06 2 . E −
06 7 . E −
07 4 .
98 5 .
01 4 .
87 4 .
90 4 .
985 10 4 . E −
05 6 . E −
05 1 . E −
04 1 . E −
04 4 . E −
05 - - - - - 615 4 . E −
06 6 . E −
06 9 . E −
06 1 . E −
05 3 . E −
06 6 .
00 5 .
83 5 .
92 5 .
94 5 . . E −
07 1 . E −
06 1 . E −
06 2 . E −
06 6 . E −
07 6 .
02 5 .
90 5 .
88 5 .
92 6 . . E −
08 9 . E −
08 1 . E −
07 2 . E −
07 5 . E −
08 6 .
01 5 .
95 5 .
82 5 .
88 6 . . E −
08 1 . E −
08 3 . E −
08 3 . E −
08 1 . E −
08 6 .
00 5 .
97 5 .
77 5 .
84 4 .
046 10 2 . E −
06 4 . E −
06 6 . E −
06 8 . E −
06 2 . E −
06 - - - - - 715 1 . E −
07 2 . E −
07 4 . E −
07 5 . E −
07 1 . E −
07 6 .
88 7 .
05 6 .
78 6 .
82 6 . . E −
08 3 . E −
08 5 . E −
08 7 . E −
08 2 . E −
08 6 .
93 7 .
04 6 .
77 6 .
82 6 . . E −
09 1 . E −
09 3 . E −
09 4 . E −
09 1 . E −
09 6 .
96 7 .
03 6 .
75 6 .
81 6 . . E −
10 2 . E −
10 5 . E −
10 6 . E −
10 1 . E −
10 6 .
99 7 .
03 6 .
75 6 .
81 6 .
997 5 3 . E −
05 4 . E −
05 7 . E −
05 9 . E −
05 3 . E −
05 - - - - - 810 1 . E −
07 1 . E −
07 3 . E −
07 4 . E −
07 1 . E −
07 7 .
88 7 .
84 7 .
80 7 .
82 7 . . E −
09 8 . E −
09 1 . E −
08 1 . E −
08 5 . E −
09 8 .
02 7 .
86 7 .
83 7 .
88 8 . . E −
10 8 . E −
10 1 . E −
09 1 . E −
09 4 . E −
10 8 .
03 7 .
92 7 .
75 7 .
82 8 . . E −
11 1 . E −
10 2 . E −
10 3 . E −
10 8 . E −
11 8 .
00 7 .
93 7 .
69 7 .
77 8 . L errors and convergence rates obtained for N ∈ { , , , , } at t end = 3 √
2. 11igure 3: 2D Lamb’s problem. Vertical velocity contour plot at t = 0 .
6. Top: ADER-DG O ×
100 elements. Bottom: reference solution computed with a P P ADER-DG schemeon an unstructured mesh of N i = 180276 triangles.in the momentum equation, (11), are chosen. We locate the source near the free surface using theDirac delta distribution at x s = (0 , −
1) and we set the related parameters to ρ s = 2200, a = − a = − ( πf c ) , f c = 14 . t D = 0 .
08. The considered homogeneous material has density ρ = 2200 andLam´e constants λ = 7509672500, µ = 7509163750, so the propagation velocities are c p = 3200 and c s = 1847 .
5. The domain is discretized using a Cartesian grid made of 200 ×
100 elements and thesimulation is run for N = 3. The vertical velocity contour plot at time t = 0 . SeisSol community code using the P P ADER-DG scheme, presented in [69, 45, 12, 13], on an unstructuredmesh made of N i = 180276 triangles. In Figure 4 we observe an almost perfect agreement between thetwo simulations in the seismogram obtained at a receiver located in x = (990 , To study the behaviour of the method under large jumps of material parameters, we consider the stiffinclusion benchmark, [75, 69]. The computational domain, Ω = [ − , × [ − . , . λ out = 2, µ out = 1 and ρ out = 1,whereas the inner material, placed in Ω in = [ − . , . × [ − . , . λ in = 200, µ in = 100 and ρ in = 1. x = (990 ,
0) obtained using ADER-DG O ×
100 elements and reference solution, computed with a P P ADER-DG scheme onan unstructured mesh of N i = 180276 triangles.The initial field is characterized by a p-wave of the form Q ( x ,
0) = r p exp (cid:20) − σ ( n · ( x − x )) (cid:21) , (29)with x = ( − . ,
0) the initial wave location, n = (1 ,
0) and, σ = 0 .
01 the standard deviation. Freesurface boundary conditions are considered, so that the normal and shear stresses vanish at the boundary.Consequently, surface waves will develop from the beginning of the simulation. At time t = 0 .
15 the planarwave reaches the stiffer material, where elastic waves propagate ten times faster that in the outer media.Then, the reflection of waves inside the inclusion yield to its vibration, which results in small amplitudewaves propagating into the outer domain. The σ xx pattern generated at time t = 0 . O ×
100 elements agrees well with the one obtained using the O ×
50 elements. To compare the results with data available in the bibliography one may referto [75] as well as to [69]. Overall, one can note a very good qualitative agreement of the wavefield withthe numerical reference solutions available in the literature.
The third test to be analysed corresponds to a 3D wave propagation problem. We consider the com-putational domain Ω = [ − , × [ − , × [ − , c p = 3200, c s = 1847 . ρ = 2200. Following [99], the initial conditionfor the vertical velocity is given by the Gaussian profile v = − . (cid:18) − r R (cid:19) , (30)13igure 5: Stiff inclusion. Component σ xx of the wavefield at time t = 0 .
3. Top: ADER-DG O ×
100 elements. Bottom ADER-DG O ×
50 elements.14igure 6: 3D wave propagation. Vertical velocity contour plot at t = 2 obtained on a mesh made of30 × ×
30 elements using the fourth order ADER-DG scheme.where R = 500 denotes the initial impulse size and r is the distance with respect to the centre of theimpulse, x = (0 , , × ×
30 elements, using the fourth order scheme. In Figure 6, the verticalvelocity contours obtained at t = 2 are plotted. To study the wave propagation three receivers havebeen placed at x r = (1000 , , x r = (3000 , − , x r = (0 , , r , the values obtained in x r for σ xx , σ xz , and u coincide with σ yy , σ yz ,and v , respectively. Similarly, at x r , σ xx matches σ yy whereas u and v take opposite values. Finally,since r is located on the x = 0 plane, we should get zero horizontal velocity. The solution obtained witha finer mesh of 60 × ×
60 elements is also included to demonstrate that mesh convergence is attained.Moreover, the obtained results are compared against a reference solution obtained with the unstructuredADER-DG code
SeisSol using a P P scheme on a mesh made of 195301 tetrahedra. Again, we cannote an excellent agreement between our solution and the numerical reference solution obtained with acommunity code. In the HSGN framework, we employ a solitary wave over a flat bottom test to assess the methodology.Next, a step-shaped bathymetry is considered and the obtained propagation of a soliton wave is comparedagainst the experimental data and the numerical results obtained for the original SGN system. Furtheranalysis on the employed ADER-DG method applied to non-hydrostatic free surface models can be foundin [5].
To assess the accuracy of the method, we study a solitary wave propagating over a flat bottom, [54, 5].The computational domain is taken to be Ω = [ − , × [ − , A = 0 . x, y ) = (0 , c = 20 and a stillwater depth H = 1. Periodic boundary conditions are considered at all boundaries. Let us remark that,for the original non-hyperbolic formulation of the SGN model, an analytical solution is available, see e.g.[14]. However, this solution does not exactly verify the hyperbolic formulation, thus it should not be15igure 7: 3D wave propagation. Seismogram at receiver x = (1000 , , O × ×
30 (red long dashed line) and M2 using 60 × × P P ADER-DG scheme onan unstructured grid of 195301 tetrahedra (black line).16igure 8: 3D wave propagation. Seismogram at receiver x r = (3000 , − , O × ×
30 (red long dashed line) and M2 using 60 × × P P ADER-DG scheme onan unstructured grid of 195301 tetrahedra (black line).17igure 9: 3D wave propagation. Seismogram at receiver x r = (0 , , O × ×
30 (red long dashed line) and M2 using 60 × × P P ADER-DG scheme onan unstructured grid of 195301 tetrahedra (black line).18mployed in a convergence study. Instead, we consider a 1D self-similar solution of the hyperbolic system(5)-(8) of the form Q ( x, t ) = Q ( ζ ) , ζ = x − V t, (31) V the velocity of the solitary wave, obtained by solving the corresponding nonlinear ODE, Q (cid:48) = ( A ( Q ) − V I ) − S ( Q ) , A ( Q ) = ∂ F ∂ Q + B ( Q ) , (32)with initial condition Q ( ζ ) = ( H , , , , (cid:15) ), (cid:15) = 10 − , H = H . The former ODE is solved using atenth order discontinuous Galerkin scheme, see [38]. The solution is used both for the initialization ofthe soliton and to compute the L errors at t end = 2. The errors and convergence rates obtained forthe density, horizontal velocity, and pressure, using polynomial degrees N ∈ { , , , , } are depicted inTable 2. Overall, the sought order of convergence is achieved but for some particular cases, in which asuboptimal order can be observed on some of the variables. Besides, Figure 10 shows the 1D profiles ofwater depth, horizontal velocity, averaged vertical velocity and pressure obtained after a one completerevolution of the soliton, t = 29 .
15. We observe that the results, computed on a mesh made of 100 × N = 5, perfectly match the initial condition. As second test case, we consider the solitary wave propagating over a step benchmark, introduced in[94]. We define the computational domain Ω = [ − , × [ − ,
1] and a step shaped obstacle of height H obs = 0 . x = 0: z b = 0 .
05 (erf(8 x ) + 1) . (33)As initial conditions, we set a soliton of amplitude A = 0 . x, y ) = ( − ,
0) and a stillwater depth H = 0 .
2. Periodic boundary conditions are imposed in x and y directions.The simulation is run using the HSGN model on a mesh made of 2000 × t ∈ { , . , . , . , . , . } . Due to the presence of the step, theamplitude of the soliton starts growing until it splits into two transmitted waves. Moreover, a reflectedwave starts to propagate in the opposite direction with respect to the soliton, followed by a train of smalldispersive waves. Notice that small spurious oscillations appear in correspondence to obstacle location x obs = 0. As already pointed out in [5], this is due to the fact that the model HSGN is rigorously validonly in the presence of a slowly varying bottom in space, which could create some problems when astrongly varying bottom topography is considered (as in the present test). To compare the numericalresults obtained with the experimental data provided in [94], we compute the ratio between the waveamplitude and the still water depth, AH = h − HH , (34)at seven different locations x ∈ {− , − , − , , , , } . In Figure 11, we observe that the results obtainedmatch pretty well the experimental data, improving the numerical results already presented in [94]. Inthe first three plots, we can observe a good agreement of the amplitude and location of the reflectedwaves, even if the already mentioned spurious oscillations at the obstacle location can be detected. Alsotransmitted waves are properly captured, including a third transmitted wave that is missing in thenumerical results in [94]. The last numerical tests aim at showing the behaviour of the proposed methodology for the simulationof the coupled problem. Two different initial water wave profiles are considered: a solitary wave and atrain of sinusoidal waves. 19 N x (cid:15) L ( h ) (cid:15) L ( u ) (cid:15) L ( p ) O L ( h ) O L ( u ) O L ( p ) Teor.3 80 1 . E −
03 8 . E −
04 8 . E −
03 - - - 4100 3 . E −
04 2 . E −
04 3 . E −
03 4 .
52 4 .
99 3 . . E −
04 9 . E −
05 1 . E −
03 4 .
81 5 .
74 4 . . E −
05 4 . E −
05 7 . E −
04 4 .
69 5 .
88 4 . . E −
05 1 . E −
05 4 . E −
04 4 .
37 5 .
86 4 .
094 60 3 . E −
04 2 . E −
04 3 . E −
03 - - - 580 6 . E −
05 4 . E −
05 7 . E −
04 5 .
72 6 .
38 5 . . E −
05 1 . E −
05 2 . E −
04 5 .
38 6 .
12 5 . . E −
06 3 . E −
06 1 . E −
04 4 .
75 5 .
89 4 . . E −
06 1 . E −
06 5 . E −
05 4 .
51 5 .
61 4 .
585 20 9 . E −
03 1 . E −
02 8 . E −
02 - - - 640 3 . E −
04 1 . E −
04 3 . E −
03 4 .
63 6 .
39 4 . . E −
05 1 . E −
05 3 . E −
04 5 .
28 5 .
51 5 . . E −
06 3 . E −
06 6 . E −
05 5 .
84 5 .
71 6 . . E −
06 9 . E −
07 1 . E −
05 5 .
89 6 .
30 5 .
726 30 3 . E −
04 3 . E −
04 3 . E −
03 - - - 740 6 . E −
05 5 . E −
05 5 . E −
04 5 .
98 6 .
25 5 . . E −
05 7 . E −
06 1 . E −
04 6 .
65 8 .
98 6 . . E −
06 1 . E −
06 2 . E −
05 9 .
17 7 .
95 8 . . E −
07 5 . E −
07 9 . E −
06 8 .
24 7 .
34 7 .
377 10 3 . E −
02 5 . E −
02 1 . E −
01 - - - 820 8 . E −
04 1 . E −
03 1 . E −
02 5 .
13 5 .
31 4 . . E −
05 5 . E −
05 6 . E −
04 6 .
07 7 .
64 6 . . E −
06 4 . E −
06 8 . E −
05 7 .
55 9 .
05 7 . . E −
07 6 . E −
07 1 . E −
05 10 .
50 8 .
45 8 . L errors and convergence rates obtained for N ∈ { , , , , } at t end = 2. 20igure 10: Solitary wave over a flat bottom. 1D profile obtained after a complete revolution of thesoliton (red dashed line) and initial solution (black line). From left top to right bottom: water depth, h ,horizontal velocity, u , averaged vertical velocity, w , pressure, p .21igure 11: Solitary wave over a step. Comparison of the time evolution of numerical results obtained usingADER-DG against the numerical and experimental data given in [94] at locations x ∈ {− , − , , , , } (from left top to right bottom). 22igure 12: Solitary wave over a step. 1D cut along y = 0 of the free surface and the bottom bathymetryat times t ∈ { , . , . , . , . , . } (from left top to right bottom).23 .3.1 Solitary wave We first study the seismic waves generated in the solid domain by the propagation of a soliton onthe water surface. We consider the computational domains Ω w = [ − , × [ − ,
2] for the fluid andΩ e = [ − , × [ − , × [ − ,
0] for the solid. As initial condition for the HSGN model we employthe planar solitary wave over a flat bottom already analysed in Section 4.2.1. The linear elasticity modelis initiated with zero values. The simulation is run using a mesh of 50 × ×
100 elements on Ω e and100 × w , (M1). Periodic boundary conditions are set in x and y directions whereas a free surfaceboundary condition is defined on the bottom of Ω e . The simulation is run up to time t end = 70, whichcorresponds to 2 . t = 58 . x r = (0 , , − x r = (40 , , − x r = (0 , , − t end = 70. The spurious oscillations generated at thebeginning of the simulation, due to the homogeneous initial condition used, quickly disappear, leading tothe smooth wave profile arising in response to the soliton. As expected, we observe that the magnitude ofthe stress tensor decreases as we get far from the sea bottom surface, while the wave front is reached at thesame time instants. On the other hand, comparison of the seismograms obtained at x r and x r provesthe conservation of the stress magnitude, as the seismic wave advances in the horizontal direction. Wealso include the results obtained with a finer grid (M2), made of 70 × ×
140 elements on Ω e and 140 × w . Finally, in Figure 15 we show the time evolution of the main variables of the non-hydrostaticmodel. Finally, we propose a sinusoidal wave propagation test. We consider the computational domains Ω w =[ − , × [ − , e = [ − , × [ − , × [ − , ∂ Q ∂t + A ( Q ) ∂ Q ∂x = S ( Q ) , A ( Q ) = ∂ F ∂ Q + B ( Q ) . (35)That is, we first decompose the conservative variables into a stationary part plus the contribution ofsmall time dependent fluctuations, i.e. Q ( x, t ) = Q ( x ) + Q (cid:48) ( x, t ). Next, we assume Q (cid:48) = ˆ Q e i ( kx − ωt ) ,where k denotes the wave number and ω is the angular frequency, obtaining the eigenvalue problem[ A ( Q ) + i E ( Q )] Q (cid:48) = ω IQ (cid:48) , E ( Q ) = ∂ S ∂ Q ( Q ) . (36)Once (36) is solved, we select one real non zero eigenvalue, λ = (cid:114) c γ + gH k + c H k + (cid:113) g H k + 2 gc H k − gc H γk + c H k + 4 c H γk + 4 c γ √ H , (37)24igure 13: Solitary wave. Solution obtained at time t = 58 .
3. Left top: 3D extrusion of the water depth, h , isosurfaces of σ zz and contour plot of σ zz at slice z = 0. From middle top to bottom right: contourplots at plane y = 0 of w , σ xz , σ xx , σ yy , and σ zz , respectively.25igure 14: Solitary wave. Seismograms at receivers x r = (0 , , −
5) (M1: red dotted line, M2: orangeline), x r = (40 , , −
5) (M1: blue dashed line, M2: cyan line), and x r = (0 , , −
20) (M1: green longdashed line; M2: light green line) obtained using ADER-DG O
4. From left top to right bottom: σ xx , σ yy , σ zz , and σ xz . 26igure 15: Solitary wave. Time evolution of the main variables of the HSGN model at pick points x p = (0 ,
0) (M1 red line) and x p = (40 ,
0) (M1 blue line) obtained using ADER-DG O
4. From left topto right bottom: h , u , w , and p . 27nd its corresponding eigenvector, and we compute the real part of the associated Q ∗ , r l = λk − H λ ( c k + gH k − λ ) tan( λt − kx )2 c k − H λ ( c k + gH k − λ )2 γc k . (38)Thus, adding the corresponding lake at rest solution, Q = ( H , , , ,
0) (39)to (38) multiplied by the sinusoidal function f ( x, y,
0) = 10 − cos (cid:18) πxs (cid:19) , (40)yields the initial condition Q HSGN ( x ,
0) = Q HSGN0 + f ( x, y, r l , (41)with H = 100 the still water depth, s = 200 the wave length, k = πs , g = 9 .
81 the gravity, and c = √ gH the celerity. Zero homogeneous initial conditions are considered for the linear elasticity model. Likewisein the previous test case, we define periodic boundary conditions on x and y directions and a free surfaceboundary condition on the bottom of Ω e . The solid domain is meshed using 108 × ×
100 hexahedralelements and the upper surface is refined with refinement factor 5 in each spatial direction to get thetwo-dimensional grid on Ω w . The simulation is run until t end = 100 using a fourth order ADER-DGscheme in both domains. The contour plots, on y = 0, of the main unknowns of the linear elasticitymodel are reported in Figure 18 for t = 50. Moreover, Figure 17 shows the 3D isosurfaces of σ zz . Theseismogram recorded at receivers x r = (0 , ,
0) and x r = (500 , , − x p = (0 , x p = (500 , In this work we have presented a high order explicit ADER-DG method for the resolution of a one-waycoupled system, describing seismic waves in the sea bottom generated by free surface water waves. Clas-sical formulations of non-hydrostatic dispersive systems and linear elasticity equations involve high ordertime and space derivatives, yielding to severe time step restrictions for explicit numerical schemes. Toovercome this issue, we propose the use of first order hyperbolic reformulations of the original systems,which lead to classical CFL restrictions with ∆ t ∝ ∆ x . In particular, we have considered the hyperbolicreformulation of the Serre-Green-Naghdi model for non-hydrostatic free surface flows and the first ordervelocity-stress formulation of linear elasticity for seismic wave propagation. Moreover, the use of hyper-bolic models allows an easier coupling of the equations and a unified discretization based on one and thesame method, i.e. employing the well known discontinuous Galerkin finite element method. High orderof accuracy in space and time has been achieved using the ADER-DG methodology based on performinga local reconstruction of the data at each cell at the aid of a space-time predictor and the subsequent28igure 16: Sinusoidal wave. Solution obtained at time t = 50. From left top to bottom right: contourplots at plane y = 0, z > − σ xx , σ yy , σ zz , σ xz , and w , respectively.29igure 17: Sinusoidal wave. Isosurfaces of σ zz , σ zz ∈ (cid:8) − , − , − , − , − , − , − , − , − , − − , − , , , , , , , , , } , at time t = 50.correction of the obtained approximation by considering the intercell flux within a classical space DGscheme. The use of a common methodology to solve both PDE systems eases the coupling between themodels, which has been done by imposing the normal stress on the upper boundary of the solid domaintaking into account the water column heigh computed using the HSGN model. The large discrepancybetween the wave lengths in the two media is addressed by considering non-conforming Cartesian gridsin the two domains. Firstly, a three-dimensional mesh for the solid domain is designed. Then, the facemesh obtained on the upper boundary of the solid domain is refined to get a mesh for the fluid do-main. A careful assessment of the developed methodology has been performed. Several benchmarks forthe linear elasticity equations and for non-hydrostatic dispersive free-surface flows are studied, showingexcellent agreement of the obtained results with available reference data. Finally, two new tests of thecoupled problem, considering solitary and sinusoidal free surface waves have been presented and allow tosuccessfully validate the proposed approach.Within this work we have assumed to have a smooth bathymetry, so that the mild bottom approx-imation, used for the derivation of the HSGN model, holds. However, practical applications might alsoinvolve non mild bottom topographies. Therefore, future research would study the use of non-hydrostaticfree-surface models for arbitrary bottom that may enlarge the applicability of the developed methodology.Besides, the simulation of non-hydrostatic flows is done using a dispersive shallow water type model, at-tending to the reduced depth of the water in comparison with the horizontal dimensions of the considereddomain. An alternative approach that might be studied, when smaller differences between the spatialdimensions are involved, is the coupling of the linear elasticity model with the fully three dimensionalfree surface Navier-Stokes equations. 30igure 18: Sinusoidal wave. Seismograms at receivers x r = (0 , ,
0) (red line) and x r = (500 , , − O
4. From left top to right bottom: σ xx , σ yy , σ zz , σ xz , u , v and w .31igure 19: Sinusoidal wave. Time evolution of the main variables of the HSGN model at pick points x p = (0 ,
0) (red line) and x p = (500 , O
4. From left top toright bottom: h , u , w , and p . 32 cknowledgements This work was financially supported by INdAM (
Istituto Nazionale di Alta Matematica , Italy) undertwo Post-doctoral grants of the research project
Progetto premiale FOE 2014-SIES ; M.D. acknowledgespartial support by the European Union’s Horizon 2020 Research and Innovation Programme under theproject
ExaHyPE , grant no. 671698 (call FETHPC-1-2014). The authors acknowledge funding fromthe Italian Ministry of Education, University and Research (MIUR) in the frame of the Departments ofExcellence Initiative 2018–2022 attributed to DICAM of the University of Trento (grant L. 232/2016) andin the frame of the PRIN 2017 project
Innovative numerical methods for evolutionary partial differentialequations and applications . Furthermore, M.D. has also received funding from the University of Trentovia the Strategic Initiative
Modeling and Simulation . References [1] P.F. Antonietti, C. Marcati, I. Mazzieri, and A. Quarteroni. High order discontinuous Galerkinmethods on simplicial elements for the elastodynamics equation.
Numerical Algorithms , 71:181–206, 2016.[2] P.F. Antonietti, I. Mazzieri, A. Quarteroni, and F. Rapetti. Non-conforming high order approx-imations of the elastodynamics equation.
Comput. Methods Appl. Mech. Eng. , 209–212:212–238,2012.[3] P.F. Antonietti, I. Mazzieri, A. Quarteroni, and F. Rapetti. High order space-time discretizationfor elastic wave propagation problems. In M. Azaiez, H. El Fekihand, and J.S. Hestaven, editors,
Proceedings of ICOSAHOM 2012, LNCSE , volume 95, pages 87–97. Springer Verlag, 2014.[4] P.F Antonietti, N. Dal Santo, I. Mazzieri, and A. Quarteroni. A high-order discontinuous Galerkinapproximation to ordinary differential equations with applications to elastodynamics.
IMA Journalof Numerical Analysis , 2017.[5] C. Bassi, L. Bonaventura, S. Busto, and M. Dumbser. A hyperbolic reformulation of the Serre-Green-Naghdi model for general bottom.
Submitted , 2020.[6] F. Bassi and S. Rebay. High-order accurate discontinuous finite element solution of the 2D Eulerequations.
J. Comput. Phys. , 138:251–285, 1997.[7] A. Bedford and D.S. Drumheller.
Elastic Wave Propagation . Wiley, Chichester, UK, 1994.[8] A. Berm´udez.
Continuum thermomechanics , volume 43 of
Progress in Mathematical Physics .Birkh¨auser Verlag, Basel, 2005.[9] A. Berm´udez, S. Busto, M. Dumbser, J.L. Ferr´ın, L. Saavedra, and M.E. V´azquez-Cend´on. Astaggered semi-implicit hybrid fv/fe projection method for weakly compressible flows.
Submitted ,2020.[10] A. Berm´udez and M. E. V´azquez Cend´on. Upwind methods for hyperbolic conservation laws withsource terms.
Computers and Fluids , 23:1049–1071, 1994.[11] J. Boussinesq. Thorie des ondes ed des remous qui se propagent le long d’un canal rectangulairehorizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareillesde la surface au fond.
Journal de Mathmatiques Pures et Appliques , 17:55–108, 1872.[12] A. Breuer, A. Heinecke, M. Bader, and C. Pelties. Accelerating SeisSol by generating vectorizedcode for sparse matrix operators.
Advances in Parallel Computing , 25:347–356, 2014.[13] A. Breuer, A. Heinecke, S. Rettenberger, M. Bader, A.A. Gabriel, and C. Pelties. Sustainedpetascale performance of seismic simulations with SeisSol on SuperMUC.
Lecture Notes in ComputerScience (LNCS) , 8488:1–18, 2014. 3314] M.O. Bristeau, A. Mangeney, J. Sainte-Marie, and N. Seguin. An energy-consistent depth-averagedEuler system: derivation and properties.
Discrete and Continuous Dynamical Systems Series B ,20:961–988, 2015.[15] L. Brugnano and V. Casulli. Iterative solution of piecewise linear systems.
SIAM Journal onScientific Computing , 30:463–472, 2007.[16] S. Busto, S. Chiocchetti, M. Dumbser, E. Gaburro, and I. Peshkov. High order ADER schemes forcontinuum mechanics.
Front. Phys. , 8:32, 2020. DOI: 10.3389/fphy.2020.00032.[17] S. Busto, J. L. Ferr´ın, E. F. Toro, and M. E. V´azquez-Cend´on. A projection hybrid high order finitevolume/finite element method for incompressible turbulent flows.
J. Comput. Phys. , 353:169–192,2018.[18] S. Busto, M. Tavelli, W. Boscheri, and M. Dumbser. Efficient high order accurate staggeredsemi-implicit discontinuous Galerkin methods for natural convection problems.
Comput. Fluids ,198:104399, 2020.[19] S. Busto, E. F. Toro, and M. E. V´azquez-Cend´on. Design and analysis of ADER–type schemes formodel advection–diffusion–reaction equations.
J. Comp. Phys. , 327:553–575, 2016.[20] M. J. Castro, J. M. Gallardo, and C. Par´es. High-order finite volume schemes based on recon-struction of states for solving hyperbolic systems with nonconservative products. applications toshallow-water systems.
Math. Comput. , 75:1103–1134, 2006.[21] V. Casulli. A semi-implicit finite difference method for non-hydrostatic free-surface flows.
Interna-tional Journal for Numerical Methods in Fluids , 30:425–440, 1999.[22] V. Casulli. A high-resolution wetting and drying algorithm for free-surface hydrodynamics.
Inter-national Journal for Numerical Methods in Fluids , 60:391–408, 2009.[23] V. Casulli. A semi-implicit numerical method for the free-surface Navier-Stokes equations.
Inter-national Journal for Numerical Methods in Fluids , 74:605–622, 2014.[24] V. Casulli and G. S. Stelling. Semi-implicit subgrid modelling of three-dimensional free-surfaceflows.
International Journal for Numerical Methods in Fluids , 67:441–449, 2011.[25] V. Casulli and P. Zanolli. Semi–implicit numerical modeling of nonhydrostatic free–surface flowsfor environmental problems.
Math. Comput. Modell. , 36:1131–1149, 2002.[26] C. Cattaneo. Sur une forme de l’quation de la chaleur liminant le paradoxe d’une propagationinstantane.
Comptes Rendues de l’Acadmie des Sciences , 247:431–433, 1958.[27] G. Chavent and B. Cockburn. The local projection p − p discontinuous Galerkin finite elementmethod for scalar conservation laws. Mathematical Modelling and Numerical Analysis , 23:565–592,1989.[28] R. Cienfuegos, E. Barthlemy, and P. Bonneton. A fourth-order compact finite volume scheme forfully nonlinear and weakly dispersive Boussinesq-type equations. Part I: Model development andanalysis.
Int. J. Numer. Methods Fluids , 51:1217–1253, 2006.[29] B. Cockburn, S. Hou, and C. W. Shu. The Runge-Kutta local projection discontinuous Galerkinfinite element method for conservation laws IV: the multidimensional case.
Math. Comput. , 54:545–581, 1990.[30] B. Cockburn, S. Y. Lin, and C. W. Shu. TVB Runge-Kutta local projection discontinuous Galerkinfinite element method for conservation laws III: one dimensional systems.
J. Comput. Phys. , 84:90–113, 1989. 3431] B. Cockburn and C. W. Shu. The Runge-Kutta local projection P1-Discontinuous Galerkin finiteelement method for scalar conservation laws.
Mathematical Modelling and Numerical Analysis ,25:337–361, 1991.[32] B. Cockburn and C. W. Shu. The local discontinuous Galerkin method for time-dependent convec-tion diffusion systems.
SIAM J. Numer. Anal. , 35:2440–2463, 1998.[33] B. Cockburn and C.W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finiteelement method for conservation laws II: general framework.
Math. Comput. , 52:411–435, 1989.[34] J. de la Puerte, M. K¨aser, M. Dumbser, and H. Igel. An arbitrary high-order discontinuous Galerkinmethod for elastic waves on unstructured meshes - IV. Anisotropy.
Geophys. J. Int. , 169:1210–1228,2007.[35] A.J.C. Barr de Saint-Venant. Thorie du mouvement non permanent des eaux, avec application auxcrues des rivires et a l’introduction de mares dans leurs lits.
Comptes Rendus de l’Acadmie desSciences , 73:147–154 237–240, 1871.[36] R. Dematt´e, V. A. Titarev, G. I. Montecinos, and E. F. Toro. ADER methods for hyperbolicequations with a time-reconstruction solver for the generalized Riemann problem: the scalar case.
Communications on Applied Mathematics and Computation , 2019.[37] F. Dhaouadi, N. Favrie, and S. Gavrilyuk. Extended Lagrangian approach for the defocusingnonlinear Schr¨odinger equation.
Studies in Applied Mathematics , pages 1–20, 2018.[38] M. Dumbser. Arbitrary high order PNPM schemes on unstructured meshes for the compressibleNavier–Stokes equations.
Computers & Fluids , 39:60–76, 2010.[39] M. Dumbser, D. S. Balsara, E. F. Toro, and C. D. Munz. A unified framework for the constructionof one-step finite-volume and discontinuous Galerkin schemes.
J. Comput. Phys. , 227:8209–8253,2008.[40] M. Dumbser, M. Castro, C. Par´es, and E.F. Toro. ADER schemes on unstructured meshes fornon-conservative hyperbolic systems: Applications to geophysical flows.
Computers and Fluids ,38:1731–1748, 2009.[41] M. Dumbser, C. Enaux, and E. F. Toro. Finite volume schemes of very high order of accuracy forstiff hyperbolic balance laws.
J. Comput. Phys. , 227:3971–4001, 2008.[42] M. Dumbser and M. Facchini. A local space-time discontinuous Galerkin method for Boussinesq–type equations.
Appl. Math. Comput. , 272:336–346, 2016.[43] M. Dumbser, F. Fambri, E. Gaburro, and A. Reinarz. On glm curl cleaning for a first order reductionof the ccz4 formulation of the einstein field equations.
J. Comput. Phys. , 404:109088, 2020.[44] M. Dumbser, F. Guercilena, S. K¨oppel, L. Rezzolla, and O. Zanotti. Conformal and covariant Z4formulation of the Einstein equations: strongly hyperbolic first–order reduction and solution withdiscontinuous Galerkin schemes.
Physical Review D , 97:084053, 2018.[45] M. Dumbser and M. K¨aser. An arbitrary high-order discontinuous Galerkin method for elastic waveson unstructured meshes II. The three-dimensional isotropic case.
Geophys. J. Int. , 167:319–336,2006.[46] M. Dumbser, M. K¨aser, and E. F. Toro. An arbitrary high-order Discontinuous Galerkin methodfor elastic waves on unstructured meshes - V. Local time stepping and p-adaptivity.
Geophys. J.Int. , 171:695–717, 2007.[47] M. Dumbser and C. D. Munz. Building blocks for arbitrary high order discontinuous Galerkinschemes.
J. Sci. Comput. , 27:215–230, 2006. 3548] M. Dumbser, I. Peshkov, E. Romenski, and O. Zanotti. High order ADER schemes for a unified firstorder hyperbolic formulation of continuum mechanics: Viscous heat-conducting fluids and elasticsolids.
J. Comput. Phys. , 314:824–862, 2016.[49] M. Dumbser, I. Peshkov, E. Romenski, and O. Zanotti. High order ADER schemes for a unified firstorder hyperbolic formulation of Newtonian continuum mechanics coupled with electro-dynamics.
J. Comput. Phys. , 348:298–342, 2017.[50] M. Dumbser, O. Zanotti, R. Loub`ere, and S. Diot. A posteriori subcell limiting of the discontinuousGalerkin finite element method for hyperbolic conservation laws.
J. Comput. Phys. , 278:47–75, 2014.[51] Michael Dumbser, Francesco Fambri, Maurizio Tavelli, Michael Bader, and Tobias Weinzierl. Effi-cient implementation of ader discontinuous Galerkin schemes for a scalable hyperbolic pde engine. axioms , 7(3):63, 2018.[52] A. Engsig-Karup, J. Hesthaven, H. Bingham, and T. Warburton. DG-FEM solution for nonlinearwave-structure interaction using Boussinesq-type equations.
Coastal Eng. , 55:197–208, 2008.[53] C. Escalante and T. Morales de Luna. A general non-hydrostatic hyperbolic formulation for Boussi-nesq dispersive shallow flows and its numerical approximation. 2020.[54] C. Escalante, M. Dumbser, and M.J. Castro. An efficient hyperbolic relaxation system for dispersivenon-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes.
J.Comput. Phys. , 394:385 – 416, 2019.[55] C. Eskilsson and S.J. Sherwin. An unstructured spectral/hp element model for enhancedBoussinesq-type equations.
Coastal Eng. , 53:947–963, 2006.[56] C. Eskilsson and S.J. Sherwin. Spectra/hp discontinuous Galerkin methods for modelling 2D Boussi-nesq equations.
J. Comput. Phys. , 212:566–589, 2006.[57] F. Fambri, M. Dumbser, S. K¨oppel, L. Rezzolla, and O. Zanotti. ADER discontinuous Galerkinschemes for general-relativistic ideal magnetohydrodynamics.
Mon. Not. R. Astron. Soc. , 477:4543–4564, 2018.[58] N. Favrie and S. Gavrilyuk. A rapid numerical method for solving Serre-Green-Naghdi equationsdescribing long free surface gravity waves.
J. Comput. Phys. , 336:104–127, 2017.[59] E.D. Fernandez-Nieto, M. Parisot, Y. Penel, and J. Sainte-Marie. A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows.
Communications in MathematicalSciences , 16:1169–1202, 2018.[60] E. Gaburro, W. Boscheri, S. Chiocchetti, C. Klingenberg, V. Springel, and M. Dumbser. High orderdirect Arbitrary-Lagrangian-Eulerian schemes on moving Voronoi meshes with topology changes.
Journal of Computational Physics , (407):109167, 2020.[61] E. Gaburro, M.J. Castro, and M. Dumbser. Well-balanced Arbitrary-Lagrangian-Eulerian finitevolume schemes on moving nonconforming meshes for the Euler equations of gas dynamics withgravity.
Mon. Not. R. Astron. Soc. , 477(2):2251–2275, 2018.[62] P. Garcia-Navarro and M.E. V´azquez-Cend´on. On numerical treatment of the source terms in theshallow water equations.
Computers & Fluids , 29:951–979, 2000.[63] G. Gassner, M. Dumbser, F. Hindenlang, and C.D. Munz. Explicit one-step time discretizationsfor discontinuous Galerkin and finite volume schemes based on local predictors.
J Comput Phys ,230:4232–4247, 2011.[64] A.E. Green, N. Laws, and P.M. Naghdi. On the theory of water waves.
Proc. R. Soc. A , 338:43–55,1974. 3665] A.E. Green and P.M. Naghdi. A derivation of equations for wave propagation in water of variabledepth.
Proc. R. Soc. A , 78:237–246, 1976.[66] M.J. Grote, A. Schneebeli, and D. Sch¨otzau. Discontinuous Galerkin finite element method for thewave equation.
SIAM J. Numer. Anal. , 44:2408–2431, 2006.[67] M. E. Gurtin.
An introduction to continuum mechanics , volume 158 of
Mathematics in Science andEngineering . Academic Press Inc., New York, 1981.[68] H. Jackson. On the eigenvalues of the ADER-WENO Galerkin predictor.
J. Comput. Phys. ,333:409–413, 2017.[69] M. K¨aser and M. Dumbser. An arbitrary high-order discontinuous Galerkin method for elasticwaves on unstructured meshes I. The two-dimensional isotropic case with external source terms.
Geophys. J. Int. , 166:855–877, 2006.[70] M. K¨aser, M. Dumbser, J. de la Puerte, and H. Igel. An arbitrary high-order Discontinuous Galerkinmethod for elastic waves on unstructured meshes - III. Viscoelastic attenuation.
Geophys. J. Int. ,168:224–242, 2007.[71] J.T. Kirby. Boussinesq Models and their Application to Coastal Processes across a Wide Range ofScales.
Journal of Waterway, Port, Coastal and Ocean Engineering , 142:1–28, 2016.[72] C. Klaij, J. J. W. Van der Vegt, and H. Van der Ven. Space-time discontinuous Galerkin methodfor the compressible Navier-Stokes equations.
J. Comput. Phys. , 217:589–611, 2006.[73] Horace Lamb. On the propagation of tremors over the surface of an elastic solid.
PhilosophicalTransactions of the Royal Society of London. Series A , 203:1–42, 1904.[74] R. J. LeVeque. Balancing source terms and flux gradients in high-resolution Godunov methods:The quasi-steady wavepropagation algorithm.
J. Comput. Phys. , 146:346–365, 1998.[75] R. J. LeVeque.
Finite Volume Methods for Hyperbolic Problems . Cambridge University Press, 2002.[76] D. Levy, C. W. Shu, and J. Yan. Local discontinuous Galerkin methods for nonlinear dispersiveequations.
J. Comput. Phys. , 196:751–772, 2004.[77] P.A. Madsen, H.B. Bingham, and H.A. Schffer. Boussinesq-type formulations for fully nonlinear andextremely dispersive water waves: Derivation and analysis.
Proceedings: Mathematical, Physicaland Engineering Sciences. , 459:1075–1104, 2003.[78] P.A. Madsen and D.R. Fuhrman.
Advances in numerical simulation of nonlinear water waves. ,chapter Higher-order Boussinesq-type modelling of nonlinear wave phenomena in deep and shallowwater. Q. Ma, Ed., World Scientific, Hackensack, NJ, 2010.[79] P.A. Madsen, R. Murray, and O.R. Sorensen. A new form of the Boussinesq equations with improvedlinear dispersion characteristics.
Coastal Eng. , 15:371–388, 1991.[80] P.A. Madsen and O.R. Sorensen. A new form of the Boussinesq equations with improved lineardispersion characteristics. Part 2. A slowly-varying bathymetry.
Coastal Eng. , 18(3):183–204, 1992.[81] G. Dal Maso, P. G. LeFloch, and F. Murat. Definition and weak stability of nonconservativeproducts.
J. Math. Pures Appl. , 74:483–548, 1995.[82] R. C. Millington, E. F. Toro, and L. A. M. Nejad.
Arbitrary High Order Methods for ConservationLaws I: The One Dimensional Scalar Case . PhD thesis, Manchester Metropolitan University,Department of Computing and Mathematics, June 1999.[83] Lucas O. Mller, Carlos Pars, and Eleuterio F. Toro. Well-balanced high-order numerical schemesfor one-dimensional blood flow in vessels with varying mechanical properties.
J. Comput. Phys. ,242:53 – 85, 2013. 3784] M.L. Munoz-Ruiz and C. Par´es. On the convergence and well-balanced property of path-conservative numerical schemes for systems of balance laws.
J. Sci. Comput. , 48:274–295, 2011.[85] O. Nwogu. Alternative form of Boussinesq equations for nearshore wave propagation.
J. Waterway,Port, Coastal, Ocean Eng. , 6:618–638, 1993.[86] C. Par´es. Numerical methods for nonconservative hyperbolic systems: a theoretical framework.
SIAM J. Numer. Anal. , 44:300–321, 2006.[87] C. Par´es and M. J. Castro. On the well-balance property of roe’s method for nonconservativehyperbolic systems. applications to shallow-water systems.
Mathematical Modelling and NumericalAnalysis , 38:821–852, 2004.[88] D. H. Peregrine. Long waves on a beach.
J. Fluid Mech. , 27(4):815–827, 1967.[89] I. Peshkov and E. Romenski. A hyperbolic model for viscous Newtonian flows.
Continuum Mech.Thermodyn. , 28:85–104, 2016.[90] W. H. Reed and T. R. Hill. Triangular mesh methods for neutron transport equation. TechnicalReport LA-UR-73-479, Los Alamos Scientific Laboratory, 1973.[91] S. Rhebergen, O. Bokhove, and J.J.W. van der Vegt. Discontinuous Galerkin finite element methodsfor hyperbolic nonconservative partial differential equations.
J. Comput. Phys. , 227:1887–1922,2008.[92] S. Rhebergen and B. Cockburn. A space-time hybridizable discontinuous Galerkin method forincompressible flows on deforming domains.
J. Comput. Phys. , 231:4185–4204, 2012.[93] S. Rhebergen, B. Cockburn, and Jaap J.W. van der Vegt. A space-time discontinuous Galerkinmethod for the incompressible Navier-Stokes equations.
J. Comput. Phys. , 233:339–358, 2013.[94] F. J. Seabra-Santos, D. P. Renouard, and A. M. Temperville. Numerical and experimental study ofthe transformation of a solitary wave over a shelf or isolated obstacle.
J. Fluid Mech. , 176:117–134,1987.[95] F. Serre. Contribution l’tude des coulements permanents et variables dans les canaux.
HouilleBlanche , 8:374–388, 1953.[96] A. H. Stroud.
Approximate Calculation of Multiple Integrals . Prentice-Hall Inc., Englewood Cliffs,New Jersey, 1971.[97] M. Tavelli and M. Dumbser. A staggered, space-time discontinuous Galerkin method for the three-dimensional incompressible Navier-Stokes equations on unstructured tetrahedral meshes.
J. Com-put. Phys. , 319:294–323, 2016.[98] M. Tavelli and M. Dumbser. A pressure-based semi-implicit space-time discontinuous Galerkinmethod on staggered unstructured meshes for the solution of the compressible Navier-Stokes equa-tions at all Mach numbers.
J. Comput. Phys. , 341:341–376, 2017.[99] M. Tavelli and M. Dumbser. Arbitrary high order accurate space-time discontinuous Galerkin finiteelement schemes on staggered unstructured meshes for linear elasticity.
J. Comput. Phys. , 366:386– 414, 2018.[100] M. Tavelli, M. Dumbser, D.E. Charrier, L. Rannabauer, T. Weinzierl, and M. Bader. A simplediffuse interface approach on adaptive Cartesian grids for the linear elastic wave equations withcomplex topography.
J. Comput. Phys. , 386:158 – 189, 2019.[101] V. A. Titarev and E. F. Toro. ADER: Arbitrary high order Godunov approach.
J. Sci. Comput. ,17:609–618, 2002. 38102] V. A. Titarev and E. F. Toro. ADER schemes for three-dimensional nonlinear hyperbolic systems.
J. Comput. Phys. , 204:715–736, 2005.[103] E. F. Toro.
Riemann Solvers and Numerical Methods for Fluid Dynamics . Springer, third edition,2009.[104] E. F. Toro, R. C. Millington, and L. A. M. Nejad.
Godunov methods , chapter Towards very highorder Godunov schemes. Springer, 2001.[105] E. F. Toro, R. C. Millington, and L. A. M. Nejad. Towards very high-order Godunov schemes.In
In Godunov Methods: Theory and Applications. Conference in Honour of S K Godunov , pages897–902, New York, Boston and London, 2001. Kluwer Academic/Plenum Publishers.[106] E. F. Toro and V. A. Titarev. Solution of the generalized Riemann problem for advection-reactionequations.
Proc. Roy. Soc. London , pages 271–281, 2002.[107] E. F. Toro and V. A. Titarev. Solution of the generalized Riemann problem for advection-reactionequations.
Proceedings of the Royal Society of London. Series A: Mathematical, Physical andEngineering Sciences , 458(2018):271–281, 2002.[108] E. F. Toro and V. A. Titarev. ADER schemes for scalar non-linear hyperbolic conservation lawswith source terms in three-space dimensions.
J. of Comp. Phys. , 202(1):196–215, 2005.[109] J. J. W. van der Vegt and H. van der Ven. Space-time discontinuous Galerkin finite element methodwith dynamic grid motion for inviscid compressible flows I. general formulation.
J. Comput. Phys. ,182:546–585, 2002.[110] H. van der Ven and J. J. W. van der Vegt. Space-time discontinuous Galerkin finite element methodwith dynamic grid motion for inviscid compressible flows II. efficient flux quadrature.
Comput.Methods Appl. Mech. Engrg. , 191:4747–4780, 2002.[111] J. Virieux. Sh-wave propagation in heterogeneous media: Velocity–stress finite–difference method.
Geophysics , 49:1933–1942, 1984.[112] Yoshiki Yamazaki, Zygmunt Kowalik, and Kwok Fai Cheung. Depth-integrated, non-hydrostaticmodel for wave breaking and run-up.
Int J Numer Methods Fluids , 61(5):473–497, 2009.[113] J. Yan and C. Shu. Local discontinuous Galerkin methods for partial differential equations withhigher order derivatives.
J. Sci. Comput. , 17:27–47, 2002.[114] J. Yan and C.W. Shu. A local discontinuous Galerkin method for KdV-type equations.
SIAM J.Numer. Anal. , 40:769–791, 2002.[115] O. Zanotti, F. Fambri, M. Dumbser, and A. Hidalgo. Space-time adaptive ADER discontinuousGalerkin finite element schemes with a posteriori sub-cell finite volume limiting.