An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on curvilinear meshes

We present an entropy stable nodal discontinuous Galerkin spectral element method (DGSEM) for the two-layer shallow water equations on two dimensional curvilinear meshes. We mimic the continuous entropy analysis on the semi-discrete level with the DGSEM constructed on Legendre-Gauss-Lobatto (LGL) nodes. The use of LGL nodes endows the collocated nodal DGSEM with the summation-by-parts property that is key in the discrete analysis. The approximation exploits an equivalent flux differencing formulation for the volume contributions, which generate an entropy conservative split-form of the governing equations. A specific combination of an entropy conservative numerical surface flux and discretization of the nonconservative terms is then applied to obtain a high-order path-conservative scheme that is entropy conservative and has the well-balanced property for discontinuous bathymetry. Dissipation is added at the interfaces to create an entropy stable approximation that satisfies the second law of thermodynamics in the discrete case. We conclude with verification of the theoretical findings through numerical tests and demonstrate results about convergence, entropy stability and well-balancedness of the scheme.


Introduction
The two-layer shallow water equations are a widely used model to describe the behavior of two immiscible fluid layers of different density, where the horizontal scales are large relative to the depth scale.Typical applications range from water bodies of different salinity and temperature [7] to mixtures of oil and water [19], as well as underwater landslides [27,37].
The system is derived from depth-averaging of the inviscid Navier-Stokes equations and can be considered as a combination of the standard shallow water equation (SWE) for each fluid layer with an additional nonconservative term that accounts for the coupling between fluid layers and interaction with the bottom topography.
The design of numerical schemes for this system encounters several difficulties, many of them connected to the inclusion of the nonconservative products.These additional coupling terms affect the wave propagation and can even change the nature of system.In contrast to standard shallow water flows, the system is only conditionally hyperbolic as the baroclinic mode allows for complex eigenvalues.This loss of hyperbolicity occurs for large velocity differences between the layers and has been linked to the onset of Kelvin-Helmholtz instability [1].The computation of the eigenvalues themselves is also not straightforward as even the onedimensional system yields a 4 × 4 flux Jacobian and a direct calculation of eigenvalues is not possible.
Another important aspect that arises in the numerical treatment of balance laws is the preservation of steadystate solutions.A varying or even discontinuous bottom topography generates additional source contributions that necessitate a careful discretization of source terms and fluxes to ensure a discrete balance.Failure to preserve these steady states can lead to accumulation of balancing errors over time and the formation of unphysical waves [9], on the order of the grid spacing.Since most shallow water applications revolve aorund the computation of small perturbations around a steady-state solution, such as gravity waves, tidal flow or tsunami formation, preserving this balance is essential.A numerical scheme that preserves these steady states is denoted as well-balanced.
In recent developments regarding numerical schemes for the two-layer SWE most activities focused on the development of finite volume methods [1,7,10,27,32,47] with only few advances in the application of high-order discontinuous Galerkin methods [15,26].Due its ability to create high-order approximations and its optimal dissipation and dispersion properties [2] discontinuous Galerkin methods can provide excellent approximations for this conditionally hyperbolic system of balance laws.In order to utilize these compelling properties, we select DGSEM as our numerical method in this paper.While the high-order nature of these methods allows for accurate approximations, the inherent low amount of numerical dissipation also makes them prone to robustness issues [21].
A contemporary approach to create robust DG approximations in the nonlinear case is the use of entropy stability [48], where the solution is bounded by a strongly convex mathematical entropy function.For the two-layer SWE such an entropy function naturally occurs as the total energy within the layers.In [21,52] a methodology to construct entropy stable DG methods based on the summation-by-parts (SBP) property and a flux-differencing formulation was proposed.Following this approach we construct an entropy stable nodal DGSEM using the SBP properties to mimic a continuous entropy analysis in the semi-discrete case.Furthemore, we show that for a special discretization of fluxes and nonconservative term the resulting scheme is well-balanced.
For the discussion in this paper we first describe the theory and challenges of the two-layer SWE in Section 2, followed by a continuous entropy analysis of the system in Section 3, which we will mimic in the semidiscrete setting.We then derive a split-form DGSEM on curvilinear elements in Section 4 that will provide a baseline scheme for the further derivation of entropy conservative and entropy stable approximations.In Section 5, we first derive entropy conservative numerical fluxes that are then used in the split-form DGSEM to create an entropy conservative approximation.We then contract the DGSEM with the entropy variables and rely on the SBP property and metric identities to demonstrate entropy conservation and well-balancedness.From this entropy conservative formulation we then extend to one that is entropy stable by adding appropriate dissipation in Section 6.The theoretical findings are then demonstrated and verified in Section 7 with a number of test cases that provide numerical results for convergence, well-balancedness and entropy conservation and stability.

Two-layer equations
We consider the two-layer SWE, which are a nonlinear system of balance laws derived from depth-averaging the inviscid Navier-Stokes-equations.The two-dimensional system is given by where the subscripts denote the upper "1" and lower "2" fluid layer.For each layer h i denotes the respective layer height, u i and v i the averaged layer velocity in xand y-direction and ρ i the different constant fluid densities, where we assume that ρ 1 < ρ 2 .Furthermore, b ≡ b(x, y) denotes the variable bottom topography and g the gravitational acceleration.The conservative part of the system corresponds to the standard SWE for each layer and only differs by additional nonconservative coupling terms on the right hand side that are given by gh 1 (h 2 ) x in the upper layer and gh 2 ( ρ 1 ρ 2 h 1 ) x in the lower layer.
To simplify the analysis we introduce a general compact notation for this system of balance laws where T denotes the state vector of conserved quantities and the nonconservative term composed of a nonconservative product and a source term is reformulated in terms of the Hadamard product with the state vector φ φ φ (u) = (0, gh 1 , gh 1 , 0, gh 2 , gh 2 ) T .To obtain a compact notation for the fluxes and nonconservative term we introduce the block vectors, see [3,52] that contain the respective components in xand y-directions The nonconservative term in (2) introduces some difficulties in the design of robust numerical methods.Additionally, in the continuous analysis, the nonconservative term renders the system only conditionally hyperbolic.This escape from the hyperbolic regime has been linked to the onset of Kevin-Helmholtz instabilities [4], where high velocity difference in the shear layer leads to heavy mixing.A mechanism that cannot be replicated in the immiscible setting.
Another difficulty arises in the computation of the eigenvalues of the flux Jacobian.Even the one dimensional system yields a fourth order characteristic polynomial, hence a direct calculation of the eigenvalues is not feasible.Instead we follow the approach from Nycander and Döös [35] and assume that the baroclinic and barotropic modes do not interact for u 1 ≈ u 2 and ρ 1 ≪ ρ 2 .This assumption leads to the eigenvalue approximation where λ ext and λ int denote the wave speeds of the barotropic and baroclinic mode, respectively written in terms of the velocities and the reduced gravity . For simplicity the onedimensional case is considered as the extension to two dimensions is straightforward.From the approximation (5) we see that internal eigenvalues may become complex and the hyperbolic regime of the system is estimated by In accordance with [17] we can get the following bound on the maximum wave speeds A particularly important property of numerical schemes for the two-layer SWE is the preservation of steadystate solutions that remain constant in time.The most important example of such a steady-state is the lake-at-rest condition As a wide range of applications are computed from small perturbations around the lake-at-rest condition, discrete preservation of ( 8) is essential.A scheme that is able to preserve ( 8) is denoted as well-balanced.

Entropy analysis
Stability estimates for numerical approximations are typically provided in terms of L 2 -stability, where the discrete solution remains bounded within the L 2 -norm, see [31,36].While L 2 -stability holds for linear problems, it fails to provide stability for general nonlinear systems [21,34].An alternative for nonlinear systems is the use of entropy stability [48], where the solution is bounded by a strongly convex mathematical entropy function.The two-layer shallow water system (1) is equipped with such a convex entropy function S = S(u), corresponding to the total energy within the layers together with the corresponding entropy flux Differentiation of the entropy function with respect to the conservative variables yields the entropy variables which are used to obtain a relation between the physical and entropy fluxes as well as the nonconservative terms Furthermore, the entropy flux potential is defined as From these compatibility relations, contracting smooth solutions of (2) with the entropy variables w yields a scalar conservation law for the entropy function We integrate the equation over the domain Ω and apply Gauss' law to rewrite the entropy flux as a surface integral and obtain the entropy conservation law (14) in integral form In the case of discontinuities entropy must be dissipated at shocks and solutions must instead satisfy the entropy inequality To prove entropy conservation of the numerical scheme, we will mimic the continuous case.Note, that the contraction of physical fluxes into entropy space depends on the chain rule, which in general does not hold discretely [40].In order to obtain entropy fluxes in the discrete setting, we therefore employ an alternative strategy based on entropy conservative numerical finite volume fluxes, that will be described in Section 5.

Split-Form DGSEM on curvinlinear elements
In this section we derive the split-form DGSEM for the nonconservative system (2), that will provide a baseline scheme from which we then develop our entropy conservative and stable approximations.To generalize the formulation to curvilinear elements we first introduce a mapping between physical and reference space.We then formulate the DGSEM for the nonconservative system and replace the volume integrals with a flux-differencing formulation to obtain the split-form DGSEM [22].

Mapping the equations
We first subdivide the physical domain Ω into K non-overlapping unstructured quadrilateral elements.Each element is mapped from physical space onto a unit square element E = [−1, 1] 2 in computational space using transfinite interpolation with linear blending [29].The mapping is given as with the coordinates in physical space → x = (x, y) T = (x 1 , x 2 ) and computational space → ξ = (ξ , η) T = (ξ 1 , ξ 2 ), respectively.From the mapping we directly calculate the covariant basis vectors that are tangential to the coordinate lines in computational space.The covariant basis vectors are then used to compute the Jacobian of the mapping as well as the contravariant basis vectors normal to the coordinate lines To map the system of equations into computational space, we first define a block matrix with the metric terms (20) together with a 6 × 6 identity matrix I to simplify the notation We transform the divergence operator into computational space using the contravariant block vectors for the flux and nonconservative terms Applying transformation (22) to (2) then yields the transformed balance law in computational space In order to obtain a free-stream preserving numerical method, the divergence of a constant flux must vanish such that the metric identities are satisfied discretely.That the metric identities hold in the discrete setting is a prerequisite for the preservation of steady-state solutions and the entropy conservation proof presented later in this work.In two space dimensions this is guaranteed for isoparametric boundaries, where the same polynomial order (or lower) from the DGSEM approximation is use to construct the mapping [28].

Nodal discontinuous Galerkin Spectral Element method
The numerical scheme developed in this work is based on a standard nodal collocation DGSEM formulation as described in the literature, see [25,29].In the following we formulate our DGSEM and extend it for nonconservative terms analogous to the numerical schemes developed in [18,42].
We begin with the derivation of a weak formulation of ( 24), where we multiply the transformed equation with test functions ϕ and integrate over the domain Ω.For nonconservative systems the definition of weak solutions is rather difficult as the theory of distributions does not apply and classical conservative approaches may recover incorrect shock speeds.Instead we follow the derivation provided by Franquet and Perrier [18], where a weak formulation is obtained from integration-by-parts for the conservative fluxes, while the nonconservative term is defined as a Borel measure according to theory of DalMaso, LeFloch and Murat [12].For smooth solutions this measure corresponds to classical integration, while at discontinuities it must be evaluated as a line integral depending on a family of paths υ.The resulting weak formulation is then given by where ŝ denotes the differential surface element given by n is the outwards pointing unit normal vector in the physical space We introduce a notation where at each interface we denote the state of the primary element as " − " and the state of the secondary (neighboring) element as " + ".Accordingly, the physical normal vector → n is defined to point from " − " to " + ", which the following relation for the normal fluxes of neighboring elements at an interface Furthermore, we introduced the inner product notation on the reference element, given for state and block vectors respectively.
To account for the definition of the nonconservative term, we further introduced the surface numerical nonconservative term ϕ ϕ ϕ T (φ φ φ • ↔ r ) ⋄ that satisfies the consistency condition and the path-conservative property [38] so that at each element interface we obtain where A NC denotes the flux Jacobian of the nonconservative subsystem and υ is a path connecting neighboring states u − and u + across a discontinuity, for complete details see [38].The numerical nonconservative term in (26) is therefore equivalent to a fluctuation form [38], where for each element (φ φ φ → n recovers the contributing fluctuations at the interface.
Next, we approximate the solution with a local tensor-product basis of polynomial degree where the basis is spanned by one-dimensional nodal Lagrange basis functions on the interval ξ ∈ [−1, 1] with N +1 interpolation nodes {ξ i } N i=0 , located at the Legendre-Gauss-Lobatto (LGL) points.In the following, interpolated quantities at the LGL nodes are denoted with capital letters and the interpolation of a function by G = I N (g).We then use the basis functions to define the discrete derivative operator apply a LGL quadrature rule to approximate the integrals present in the weak formulation (26) and collocate the interpolation and quadrature nodes.This yields the diagonal mass matrix with the LGL quadrature weights {ω i } N i=0 .Using the Kronecker delta property of the Lagrange basis functions (34) we then introduce the discrete inner product notation This specific choice of interpolation and differentiation operators possesses the diagonal norm SBP property for all polynomial orders N [20] where we introduced the SBP matrix Q and the boundary matrix This property will be necessary to prove entropy stability, as it allows us to mimic integration-by-parts and enables the extended Gauss Law in the discrete case [30].We introduce polynomial approximations for the solution variables, fluxes, test functions and the nonconservative terms as well as LGL quadrature to obtain a discrete version of (26) where capital letters denote the discrete approximation and the surface fluxes We then mimic integration-by-parts with the the SBP property (38) and apply the discrete extended Gauss Law from [30] on the flux terms to obtain the strong form DGSEM approximation of ( 24) The method (42) serves as a baseline from which we will continue to construct an entropy stable approximation using a special combination of differencing operator, numerical flux, and numerical nonconservative term.

Flux-differencing formulation
To obtain an entropy stable discretization, we make use of another property that diagonal norm SBP operators have and rewrite the flux divergence term in the volume into a flux differencing formulation [5,22] with the consistent and symmetric two-point flux F # and arithmetic mean values given by This formulation was originally found by Fisher et al. [16], who showed that a diagonal norm SBP operator can be rewritten into a finite volume type flux-differencing formulation on a staggered grid.From a particular choice of the two-point volume flux, the formulation (43) recovers different alternative split-forms of the original partial differential equation that can improve robustness of the scheme [22].The formulation in (43) also extends two-point finite-volume fluxes to high-order [22] and recovers an entropy conservative volume discretization if an entropy conservative two-point flux is used [5,21].Recently, this approach has been applied successfully to create entropy stable DG methods for the shallow water equation [23,50,53], magnetohydrodynamics equations [3,44] and the compressible Navier-Stokes equations [8,52].
In [42,49] it was shown that the flux-differencing formulation can also be extended to nonconservative systems such that we can apply the same discretization technique on the conservative and nonconservative parts of the equation.This ensures that both terms are evaluated in the same way, which will be necessary to obtain entropy conservation and discrete preservation of steady-state solutions.Therefore, we introduce the following flux-differencing formulation for the nonconservative volume term where (Φ Φ Φ • ↔ R) # denotes a two-point volume numerical nonconservative term with a given structure where Φ Φ Φ # and ↔ R # are vector-valued functions.It was shown in [49,Lemma 2] that under the consistency conditions the volume numerical nonconservative term (46) provides a consistent volume discretization of the nonconservative term.For more details about the flux-differencing formulation for nonconservative terms we refer to [42,49].
Inserting the flux-differencing formulations ( 43) and ( 45) into the discrete strong formulation (42) then yields the split-form DGSEM With method (48) we are now equipped to derive the necessary components ⋄ to first create an entropy conservative and well-balanced scheme, and from that, a scheme that is entropy stable.

Entropy-Conservative DGSEM
To show that the scheme is entropy conservative we use the approach described in [3,50,52] and mimic the continuous entropy analysis of Section 3 in the discrete setting.First, we contract the split-form DGSEM (48) into entropy space by setting the test function to the interpolated entropy variables W

Time derivative
As we are interested in the semi-discrete analysis, we assume that the chain rule with respect to differentiation in time holds and examine the time derivative term in (49) to obtain From (50) we see that the rate of entropy change for a single element depends solely on the remaining surface and volume contributions.To obtain the total entropy, we sum over the all K elements in the domain To demonstrate discrete entropy conservation, we will show that the entropy conservation law ( 14) is satisfied in a discrete sense by the split-form approximation (48).Assuming a closed system, the second law of thermodynamics states that the entropy flux at the boundaries vanishes and the total entropy must be conserved

Discrete entropy flux
The contraction of the physical flux nonconservative terms to the entropy flux in ( 12) is challenging as it depends on the chain rule, which in general does not hold in the discrete setting [40].Therefore, we apply an alternative strategy to find a combination of numerical fluxes and numerical nonconservative terms that recover the entropy flux on a discrete level.To reduce the degrees of freedom in this derivation we assume a fixed discretization of the numerical nonconservative term which can be derived from the path-conservative property, using a linear path and satisfies the consistency conditions (47).The task now reduces to finding an accompanying entropy conservative (EC) numerical flux To derive such EC fluxes we require discrete entropy conservation in a finite volume context similar to [13,51].These finite volume type numerical fluxes are then directly applicable to the split-form DGSEM (43) as the volume contribution relates to a subcell finite volume scheme and elements are coupled via numerical fluxes on the surface [16,22].We consider a single interface between neighboring finite volume cells with distinct cell averages denoted by " − " and " + " and cell size ∆x.We introduce a notation to define jumps and arithmetic averages of cell values across the interface The respective finite volume formulation of system (2) with the numerical nonconservative term ( 53) is given by at each side of the interface.To obtain an entropy conservation statement we follow the continuous analysis and contract (55) with the entropy variables.As before, we assume time continuity such that the chain rule in time holds to obtain the rate of entropy change in each cell We then sum over both cells to find the total discrete entropy change over the interface From (57) we see that a discrete version of the integral entropy conservation law (15) holds provided the right hand side recovers the jump of the entropy flux over the interface Manipulating (58) and subsequently substituting the definition of the entropy potential (13) we find an entropy conservation condition for the numerical flux Accordingly, a numerical flux that satisfies this condition in combination with the numerical nonconservative term (53) recovers the integral entropy conservation law (15) and is therefore defined as an EC numerical flux.
With this in mind we create an EC split-form DGSEM by selecting a proper combination of EC fluxes and nonconservative terms for both volume and surface contributions.In general the split-form offers flexibility in choosing different discretizations on the volume and the surface, which was used in [50] to obtain a wellbalanced and EC scheme for the standard SWE.
However, we find that with the numerical nonconservative term (53) both properties can be achieved using the same numerical nonconservative term and EC flux for the volume and surface parts in (49).The numerical nonconservative term (53) can also be applied in the volume as it corresponds to the structure given in (46) when setting and satisfies the consistency conditions (47).This circumvents the need for different surface and volume fluxes, such that our method is given by We then set the following EC numerical fluxes which fulfill the entropy conservation condition (59) as shown in Appendix A for F EC 1 .This EC flux (62) corresponds to the EC volume flux found by Wintermeyer et al. [50] for the standard SWE within each layer.The reason behind this is that the conservative part of the twolayer model (1) corresponds to the standard shallow water system in each layer and only the nonconservative terms differ.Remark 5.1.We note that due to this similarity the EC split-form DGSEM (61) can be directly applied to the standard SWE if the EC flux from [50] and the numerical nonconservative term (53) are used.
To demonstrate entropy conservation for approximation (61) we next examine the volume and surface contributions separately.We start with the volume parts and show that for the given choice of two-point volume flux (62) and numerical nonconservative term (53) they become the entropy flux evaluated at the boundary, when contracted into entropy space.Using the same combination of EC flux and numerical nonconservative term the total surface contribution then cancels in entropy space, which shows entropy conservation.Furthermore, we demonstrate that by using an analogous discretization for the pressure and nonconservative terms yields a method that is well-balanced.

Volume contribution
Lemma 5.2 (Entropy contribution of the curvilinear volume terms).The curvilinear volume contributions of the split-form DGSEM (48) with the EC volume flux (62) contracted into entropy space become the entropy flux evaluated at the interface.
Proof.See Appendix B.
The result from Lemma 5.2 shows that the volume contributions from the flux and nonconservative terms cancel in entropy space and become the entropy flux evaluated at the boundary.Therefore, the volume parts do not contribute to the total discrete entropy and the entropy balance is entirely dependent on the surface contributions.

Surface contribution
On the surface we select the same EC numerical flux as in the volume.From Lemma 5.2 we know that the volume contributions move onto the surface such that the change in discrete entropy depends solely on the surface terms.So, we assemble the surface contributions for a single element To obtain the total discrete entropy as described by (51) we sum over all elements k = 1, ..., K. Due to the discontinuous approximation, summing over the elements creates jumps in the fluxes, nonconservative terms and entropy variables, whereas the numerical surface fluxes are defined uniquely at each interface.
After summing over the elements, we examine the terms in the surface integrals (64) separately.Starting with the physical flux term we obtain jumps in the flux and entropy variables and use that the numerical flux is unique at the interface to have In a similar way summing over the second integral generates a jump in the entropy flux The third integral in (64) contains the the numerical nonconservative term (53), which is non-unique at the interface.In the following we examine the total surface contribution of the nonconservative parts after summing over the interfaces.
Lemma 5.4 (Surface contribution of the numerical nonconservative term).Setting the numerical nonconservative term to (53) yields the total surface contribution of the nonconservative terms Proof.The total surface contribution from the nonconservative terms in (64) is given by Introducing the definition for the numerical nonconservative term ( 53) into (64) and examining the nonconservative contribution, where " − " denotes the primary and " + " the secondary state, yields To obtain the total interface contribution, we must take contributions from neighboring cells into account.Since the nonconservative term is not symmetric, we need to consider the different normal directions of neighboring cells as defined in (29).Summing over all elements, we then obtain the total contribution from the interfaces where we used that due to anti-symmetry the different sign in the jump operator between " + " and " − " cancels with the sign flip from opposite normal directions.
Lemma 5.5 (Total entropy contribution of the surface terms).When summing over all elements, the total entropy contributions of the advective and nonconservative terms in (64) cancel on the surface Proof.First we write the fluxes in terms of jumps according to (65), (66) and apply the result of Lemma 5.4 for the nonconservative contributions to obtain Now we substitute the definition of the entropy flux potential ( 13) and substitute the entropy conservation condition (59) to cancel the remaining terms and therefore show that the total entropy contribution on the surface is zero.Now that the surface and volume contributions are examined, we have everything assembled and combine the previous results to show that approximation (48) is EC.
Theorem 5.6 (Entropy conservation of the curvilinear split-form DGSEM for the two-layer shallow water equations).The curvilinear split-form DGSEM for the two-layer shallow water equations with the EC flux (62) and the numerical nonconservative term (53) is entropy conservative.
Proof.From Lemma 5.2, we have that the entropy contributions cancel in the volume and generate the entropy flux on the surface.Then we apply the result from Lemma 5.5 to show that the entropy contribution cancels on the surface and therefore entropy is conserved.

Well-Balancedness
In the following we show that in addition to entropy conservation, the discretization ( 74) is well-balanced as it satisfies the lake-at-rest condition (8).
Corollary 5.7 (Well-balancedness of the curvilinear split-form DGSEM for the two-layer shallow water equations).The EC curvilinear split-form DGSEM for the two-layer shallow water equations preserves the lake-at-rest initial condition (8) Proof.To show that the scheme is well-balanced, we take the same approach as demonstrating entropy conservation.That is we examine the volume and surface parts separately, to see that both contributions vanish for the lake-at-rest condition (8).
In the continuity equations in (4), it is clear that well-balancedness is directly satisfied as all products vanish due to the initial condition u 1 , u 2 , v 1 , v 2 = 0. Therefore, only the momentum equations need to be considered to show well-balancedness.For simplicity the proof is only demonstrated for the h 1 u 1 -equation as the other components are analogous.
Inserting the lake-at-rest conditions (8) into the EC approximation (74) and using that the volume flux contribution simplifies to The nonconservative volume contribution is then rewritten into a similar form.Therefore, we first use that the metric identities ( 25) are satisfied such that any terms containing only local contributions By expanding the nonconservative term, the contribution to the h 1 u 1 -equation is then given by Next we assemble the total volume contribution from (77) and (79).We first use that b + h 1 + h 2 = constant, according to (8).Then, from the consistency of the derivative matrix and assuming that the metric identities (25) hold discretely, the volume parts vanish for the lake-at-rest condition The remaining surface contributions are given by Expanding the terms in the h 1 u 1 -equation, again all products with velocity components vanish, and (81) simplifies to The first three terms are rewritten as Inserting this into (82) and using that b + h 1 + h 2 = 0 according to the lake-at-rest condition yields which shows that the h 1 u 1 -equation satisfies (8) as both the volume and the surface contributions vanish.As the remaining momentum equations are analogous, this shows that the EC split-form DGSEM approximation (74) is well-balanced for the two-layer SWE.

Entropy-Stable DGSEM
The high-order approximation introduced in Theorem 5.6 satisfies the entropy conservation statement ( 14) in a semi-discrete setting.In the presence of discontinuities entropy conservation is not sufficient and instead the entropy inequality ( 16) must be satisfied to account for entropy dissipation at shocks.Following the idea of Tadmor [48], we add numerical dissipation to the EC scheme from Theorem 5.6 to construct a discretization that is entropy stable (ES).We use the approach from [17] and create an ES flux by adding a local Lax-Friedrichs type dissipation to the EC flux where λ max denotes an approximation of the maximum eigenvalue from ( 7) and H = w u is the symmetric positive definite entropy Jacobian matrix, evaluated at the arithmetic average value of the solution state at element interfaces.That the dissipation is dependent on the jump in entropy variables guarantees that entropy is dissipated at shocks and conserved for smooth solutions.
Theorem 6.1 (Entropy stability of the curvilinear split-form DGSEM for the two-layer shallow water equations).
The curvilinear split-form DGSEM for the two-layer shallow water equations is ES and well-balanced, when the ES flux (85) is used at element interfaces, together with the EC flux (75) on the volume and the numerical nonconservative term (53).
Proof.Entropy stability follows, provided (86) is more dissipative than the EC discretization (74).From Lemma 5.5 it is clear that the surface contributions of (86) vanish together with the EC flux part of (85) when contracted into entropy space.Therefore, the change in total entropy obtained from summing over all elements depends only on the remaining dissipation term Using that H is the Hessian of a convex function and therefore positive definite, the total entropy decreases in time and the discretization is therefore ES.For the lake-at-rest initial condition (8), the entropy variables (11) remain constant and the ES flux simplifies to the EC flux.The proof for well-balancedness then follows directly from Corollary 5.7.

Results
To demonstrate the theoretical findings and performance of the split-form DGSEM, we present and discuss results from several numerical experiments.First, we investigate convergence properties and show spectral convergence on a curvilinear mesh.Then we provide numerical evidence of the well-balanced property for discontinuous bottom topography and add a perturbation to this test case to demonstrate entropy stability.We then close the discussion with a more complex test, showcasing a partial dam break of a parabolic dam.Numerical results are obtained using the open-source framework Trixi.jl[41,46].A reproducibility repository is available on Zenodo [14] and GitHub1 for the results presented herein.All computations use a low storage five-stage fourth-order Runge-Kutta scheme of Carpenter and Kennedy [6] for time integration, with either fixed or CFL-based time stepping [11].

Convergence Results
We apply the method of manufactured solutions to demonstrate spectral convergence of the numerical approximation against an exact reference solution.To construct the exact solution, layer heights and bathymetry are defined by trigonometric functions and constants are used for the velocities u 1 = 0.9, u 2 = 1, v 1 = 1, v 2 = 0.9 b = 1 + 0.1 cos(πx) sin(πy) The additional source terms for the manufactured solution are then computed using analytical derivatives with the gravitational constant set to g = 10 and the density ratio ρ 1 ρ 2 = 0.9.Spectral convergence for both the ES and EC approximations is shown in Figure 2, which shows a semi-log plot for the error in L 2 over the polynomial degree N. We see suboptimal convergence for the EC scheme for various polynomial degrees N. The characteristic pattern of suboptimal convergence at odd polynomial degree, reported in [20,22] cannot be observed for curved elements, but does occur under the initial conditions (88) and resolution parameters on a mesh with straight element edges.The convergence behavior of the remaining quantities was similar and is not shown.

Well-balanced
The next test is adapted from [50] and demonstrates that the approximation preserves the lake-at-rest condition (8) for the two-layer SWE with discontinuous bottom topography, according to Corollary 5.7.Again we use the curvilinear mesh from Figure 1 and set the following lake-at-rest conditions with a discontinuous bottom topography at element 3 × 2 (see Figure 3) We set periodic boundary conditions and compute the solution with a polynomial degree N = 8 and a CFLbased timestep up to the final time t end = 100.In Figure 3 we show the bottom topography and present the error in the absolute water height H 1 = h 1 + h 2 + b relative to the lake-at-rest initial condition (90).Even after considerate time of t = 100, errors in H 1 and H 2 remain around machine precision and the steady-state is preserved.Similar results with errors near unit roundoff were also obtained with the ES flux.

Well-balanced Perturbation
Next, we introduce a small discontinuous perturbation to the lake-at-rest configuration described in Section 7.2, and examine the resulting change in entropy for both the EC and ES schemes.We adopt the previous setup, but To demonstrate entropy conservation for the approximation (74) and entropy stability for (86), we evaluate the total discrete entropy change within the domain Ω as Minimum, mean and maximum values of the entropy change for ↔ F EC and ↔ F ES within the time interval T = [0, 0.1] are presented in Table 1.For the EC flux we see that entropy is conserved discretely with a total discrete entropy change around machine precision.For the ES flux on the other hand, we observe a noticeable change as entropy is dissipated at discontinuities.The solutions shows a strict decrease in entropy throughout the time interval T and therefore satisfies the entropy inequality (16).

Parabolic dam break
As a final test case, we apply the EC and ES schemes to a more complex problem with non-periodic BCs and demonstrate the solution behavior.We consider a partial dam break problem with initial conditions similar to [33,39], but with a modified geometry featuring a parabolic dam to create a curvilinear boundary.The solution domain is given by Ω = [0, 10]2 , and includes a parabolic dam described by the centerline x = 1 25 y 2 − 0.4y + 6, thickness d = 0.2 and a gap size of 1.We solve the problem on a curvilinear quadrilateral mesh containing 2239 elements and set the boundary conditions as slip-wall for the dam and Dirichlet for the outer boundaries.The curvilinear quadrilateral mesh was constructed using HOHQMesh.jl 2 .Initial conditions are chosen such that both layers are at rest, with a discontinuity across the dam.
In Figure 4 a sequence of computational results for the ES scheme obtained with polynomial degree N = 3 is shown over time.After the initial dam break we observe a shock and rarefaction wave in both layers and the formation of vortices at the corners of the dam opening.Furthermore, in Figure 5 we compare the result with the EC scheme without dissipation along the horizontal line y = 5.We see that while both schemes do well in resolving the shock and rarefaction waves.The dissipation-free EC discretization shows significant spurious oscillations, due to missing entropy dissipation at shocks.The ES scheme on the other hand dissipates at shocks and eliminates notable oscillations.However, the additional dissipation is only designed to satisfy an entropy inequality and additional shock capturing is necessary to obtain a scheme that is oscillation free, while remaining ES, eg.[24,43] 4.00 4.

Conclusion
We presented a high-order collocated nodal discontinuous Galerkin spectral element method (DGSEM) for the two-layer shallow water equations on curvilinear quadrilateral meshes.We first constructed the DGSEM in a specific way such that it is endowed with the summation-by-parts (SBP) property and introduced a path-conservative approximation of the nonconservative terms.Using the SBP property we then replaced the volume contributions with a flux-differencing formulation and demonstrated that the resulting approximation is entropy conservative (EC) for a specific combination of the EC flux and numerical nonconservative terms.Applying an equivalent discretization for the pressure and nonconservative terms, we further show that the method is wellbalanced for discontinuous bottom topography.From the dissipation free EC formulation we then constructed an entropy stable scheme by adding additional numerical dissipation at interfaces using a local Lax-Friedrichs type dissipation as a closed form of the eigenvalues is not available.Finally, we provide numerical results for a number of academic test cases to demonstrate convergence, wellbalancedness and entropy stability of the scheme.We then conclude with the application on a more complex parabolic dam break test and demonstrate the solution behavior.The numerical tests verify the analysis and show that the scheme is well-balanced and EC up to machine precision.
In future work, we aim to complement the scheme with wetting and drying techniques and add shock capturing methods to create a scheme that is oscillation-free.Furthermore, we aim to extend the present formulation to obtain an entropy stable and well-balanced discretization of the Savage-Hutter model [45] for submarine avalanches.
is analogous.The contraction of the EC flux in the ξ -direction is given by

B Entropy Volume Contribution
To show that the volume contributions generate entropy fluxes at the boundary, we will examine the volume contributions given by First we expand the volume contributions of the fluxes and the nonconservative terms We will only present the proof for the ξ -direction as it simplifies the analysis and the η-direction is done analogously.Next, we rewrite (101) in terms of the undivided differencing operator (103) We make use of the SBP property 2Q im = Q im − Q mi + B im and subsequently exchange the indexing of i and m to rewrite Q mi as Q im , using that ↔ F EC (U i j , U m j ) and J → a 1 (i,m) j are symmetric regarding i and j We then use consistency with the physical flux

n
and numerical nonconservative terms (Φ Φ Φ • ↔ R) ⋄ • → n are replaced by the surface normal fluxes F n and surface numerical nonconservative terms (Φ Φ Φ • ↔ R) ⋄ n .Since the physical fluxes are not uniquely defined at the discontinuous interfaces, we replace the surface normal flux F n by a surface numerical flux F * n , to resolve the discontinuity such that we obtain the discrete weak formulation

Figure 1 :
Figure 1: Curvilinear mesh adapted from [50] with polynomial degree N = 6 representation of interfaces used to test spectral convergence and well-balancedness The resulting problem is solved in the domain [0, √ 2] 2 on the curvilinear mesh shown in Figure 1 using polynomial degree N = 6 representation at element edges.Solutions are computed at the final time t end = 1 for polynomial degrees up to N = 30 and a fixed timestep ∆t = 1/12000.Spectral convergence for both the ES and EC approximations is shown in Figure2, which shows a semi-log plot for the error in L 2 over the polynomial degree N. We see suboptimal convergence for the EC scheme for various polynomial degrees N. The characteristic pattern of suboptimal convergence at odd polynomial degree, reported in[20,22] cannot be observed for curved elements, but does occur under the initial conditions (88) and resolution parameters on a mesh with straight element edges.The convergence behavior of the remaining quantities was similar and is not shown.

Figure 2 :
Figure 2: Spectral convergence of ES and EC approximations in space shown in a semi-log plot for the L 2 -Error in h 1 u 1 over the polynomial degree N. Results are obtained with a fixed timestep of ∆t = 1/12000 at the final time t end = 1

Figure 3 :
Figure 3: (a) Discontinuous bottom topography in element 3 × 2 and (b) lake-at-rest error at time t end = 100 for the EC flux with N = 8 and CFL = 0.7

Figure 5 :
Figure 5: Layer heights along y = 5 for the EC and ES flux at t = 0.25 with polynomial degree N = 3 and ∆t = 10 −4

F 1 (F
EC (U i j , U m j ) • J → a EC (U i j , U m j ) • J → a 1 (i,m) j .

↔ F EC = ↔F↔F
at the boundary part, as the boundary matrix consists of zero entries besides B 00 = −1 and B NN = 1 to substitute the definition of the entropy flux potential (13)B im W T i j EC (U i j , U m j ) = B im → Ψ i j + →