On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows

In this paper we propose a new reformulation of the first order hyperbolic model for unsteady turbulent shallow water flows recently proposed in Gavrilyuk et al. (J Comput Phys 366:252–280, 2018). The novelty of the formulation forwarded here is the use of a new evolution variable that guarantees the trace of the discrete Reynolds stress tensor to be always non-negative. The mathematical model is particularly challenging because one important subset of evolution equations is nonconservative and the nonconservative products also act across genuinely nonlinear fields. Therefore, in this paper we first consider a thermodynamically compatible viscous extension of the model that is necessary to define a proper vanishing viscosity limit of the inviscid model and that is absolutely fundamental for the subsequent construction of a thermodynamically compatible numerical scheme. We then introduce two different, but related, families of numerical methods for its solution. The first scheme is a provably thermodynamically compatible semi-discrete finite volume scheme that makes direct use of the Godunov form of the equations and can therefore be called a discrete Godunov formalism. The new method mimics the underlying continuous viscous system exactly at the semi-discrete level and is thus consistent with the conservation of total energy, with the entropy inequality and with the vanishing viscosity limit of the model. The second scheme is a general purpose high order path-conservative ADER discontinuous Galerkin finite element method with a posteriori subcell finite volume limiter that can be applied to the inviscid as well as to the viscous form of the model. Both schemes have in common that they make use of path integrals to define the jump terms at the element interfaces. The different numerical methods are applied to the inviscid system and are compared with each other and with the scheme proposed in Gavrilyuk et al. (2018) on the example of three Riemann problems. Moreover, we make the comparison with a fully resolved solution of the underlying viscous system with small viscosity parameter (vanishing viscosity limit). In all cases an excellent agreement between the different schemes is achieved. We furthermore show numerical convergence rates of ADER-DG schemes up to sixth order in space and time and also present two challenging test problems for the model where we also compare with available experimental data.


Introduction
In the last decades, a lot of work has been devoted to the study of shallow water flows. When dispersive effects are negligible (this is the case, for example, for the modelling of hydraulic jumps for large Froude numbers, or tsunami waves), one usually employs the classical Saint-Venant (SV) or shallow water equations. The underlying hypothesis in the derivation of the Saint-Venant equations is the flow potentiality. The horizontal vorticity (parallel to the bottom) in the shallow water approximation is related with the horizontal velocity shear: ω || ≈ V z , where V is the instantaneous (non-averaged) horizontal velocity, and the index z means the derivative in the vertical direction. The absence of the vorticity means the absence of the horizontal velocity shear. The shallow water equations are hyperbolic, see e.g. [116]. When discontinuous solutions to the SV equations are studied, one uses the conservation of mass and momentum at the shocks, but the energy equation is always used as the entropy inequality. The reason for this is the following. The fluid flow ahead of the jump front is supercritical with respect to the front and almost potential, while behind the front it is highly turbulent: large vortex structures are usually formed. The energy is not conserved because a part of this energy is transformed into the energy of vortexes which is not taken into account in the SV model. An ideal model of free surface shallow flows which takes into account shear effects was recently derived by Teshukov [112]. The governing equations are obtained by depth averaging of the multi-dimensional Euler equations [100,101,112]. The hypothesis of smallness of the horizontal vorticity (the hypothesis of weakly sheared flows) allows us to keep the second order depth averaged correlations in the governing equations but neglect the third order correlations, and thus to close the governing system in the dissipationless limit. To apply the model to the study of real flows (formation of roll waves and hydraulic jumps) the model was complemented by dissipative source terms, see [69,80,100,101].
The corresponding multi-dimensional model of shear shallow water flows is a hyperbolic system of equations which is reminiscent of the Reynolds-averaged Euler equations for barotropic compressible turbulent flows. The model has three families of characteristics corresponding to the propagation of surface waves, shear waves and waves propagating with the average flow velocity. The main difficulty in studying such a system is the highly non-conservative nature of the governing equations: for six unknowns (the fluid depth, two components of the depth averaged horizontal velocity, and three independent components of the symmetric Reynolds stress tensor) one has only five conservation laws: conservation of mass, momentum, energy and mathematical "entropy". The last one determines the evolution of the determinant of the Reynolds stress tensor. The non-conservative nature of the multi-dimensional equations represents an enormous difficulty from the mathematical and numerical point of view. The definition and computation of discontinuous solutions for non-conservative hyperbolic equations is a challenging problem, see e.g. [5,23,85]. A numerical method (based on a splitting procedure) was recently developed for solving this non-conservative system [69,80]. The essential ingredient was the use of the energy conservation. It allowed, in particular, creation of vorticity once the jump appears. The splitting procedure was as follows. First, a geometric splitting was applied consisting in solving the governing equations first in x and then in y direction. Second, each one-dimensional system was also split into two subsystems, each of which contained only one 'sound' speed: the velocity of surface waves for the first sub-system, and the velocity of shear waves for the second sub-system. Each subsystem admitted its own energy conservation law, and its own "entropy". However, such an operator splitting could be also a source of numerical errors. This is why it is very important to develop also different numerical methods for solving this challenging non-conservative system, like the two new unsplit schemes proposed in this paper, namely a completely new thermodynamically compatible unsplit finite volume scheme, as well as a slightly modified general-purpose high order ADER discontinuous Galerkin finite element method.
In order to put the new unsplit thermodynamically compatible finite volume scheme presented in this paper into the proper context, let us briefly review the main ideas on which it is based. Exactly sixty years ago, in 1961, Godunov published his groundbreaking paper An interesting class of quasilinear systems [70], in which he discovered the connection between symmetric hyperbolicity in the sense of Friedrichs [62] and thermodynamic compatibility (SHTC), ten years before Friedrichs and Lax [63], who independently rediscovered the same connection again in 1971. In a subsequent series of papers, Godunov and Romenski carried out further research on this link between symmetric hyperbolicity and thermodynamic compatibility and generalized the seminal idea of Godunov to the more general SHTC framework of symmetric hyperbolic and thermodynamically compatible systems, which includes not only the compressible Euler equations of gasdynamics, but also the magnetohydrodynamics (MHD) equations [71] and the equations of nonlinear hyperelasticity [74,75,77]. The findings of Godunov and Romenski on nonlinear hyperelasticity were subsequently further employed and extended in [1,6,14,17,47,57,59,67,72,81,87,88,91,93]. A very general class of symmetric hyperbolic and thermodynamically compatible systems was presented by Romenski in [102], which is able to describe the interaction of moving-dielectric solids with electromagnetic fields, the dynamics of superfluid helium and also contains a hyperbolic model for heat conduction. An extension of this class of models to compressible multi-phase flows was forwarded in [103,104,106]. The SHTC framework remains valid even in the context of special and general relativity, see [76,105]. Recently, a connection between the class of symmetric hyperbolic and thermodynamically compatible systems and Hamiltonian mechanics was rigorously established in [92]. SHTC systems go also beyond classical continuum mechanics, see e.g. [94] for an SHTC formulation of continuum mechanics with torsion. Despite the mathematical beauty and rigor of the SHTC framework, up to now it was never carried over to the discrete level. So far, most papers on on thermodynamically compatible schemes are based on the ideas of the seminal work of Tadmor [108], in which a discrete extra conservation law for the entropy is obtained as a consequence of the discretization of all other equations (including the energy conservation, which is explicitly discretized). Instead, in the new scheme presented in this paper we are not discretizing the energy equation explicitly, but are rather looking for a thermodynamically compatible scheme in which a discrete total energy conservation law is obtained as direct consequence of the compatible discretization of all the other equations. For an interesting application of entropy compatible schemes for the discretization of non-conservative equations, see [3,60]. The ideas presented there are related to the new compatible scheme introduced in the present paper, though [3,60] deal with much simpler equation systems. Recently, convergence of entropy-stable schemes was proven in [26]. For extensions of entropy-compatible schemes to high order discontinuous Galerkin methods, see [28,36,66,78,84] and references therein. While most of the aforementioned schemes are thermodynamically compatible only at the semi-discrete level, a fully discrete entropy-stable scheme has been recently presented in [95]. We also would like to point out that a very general framework for the construction of numerical schemes satisfying additional extra conservation laws has been recently forwarded by Abgrall in [2].
As already stated above, the major difference of the thermodynamically compatible scheme proposed in this paper with respect to previous thermodynamically compatible schemes is its discrete compatibility with the conservation of total energy as a consequence of all equations and not the conservation of entropy as a consequence. In other words, the thermodynamically compatible finite volume scheme presented in this paper never explicitly discretizes the energy equation, but total energy conservation is obtained as a mere consequence of a thermodynamically compatible discretization of the other equations, including a compatible discretization of the numerical viscosity.
The second unsplit scheme proposed in this paper is a fully-discrete one-step high order ADER discontinuous Galerkin method (ADER-DG). Explicit discontinuous Galerkin schemes for hyperbolic equations have been put forward by Reed and Hill in [96] introducing the use of discontinuous polynomials in a Galerkin framework to allow the jump of the discrete solution across cell boundaries.
Then, the first extensions to multidimensional and non linear hyperbolic systems were presented in the series of papers by Cockburn and Shu [27,[30][31][32][33]. Parabolic terms have been considered for the first time in [9,10,34,35]. The severe time step restriction induced by the inclusion of higher order derivatives, [83,122,123], and nonlinear dispersive equations, [51,53,54], has driven to the development of fully implicit approaches, [44], whose major disadvantage is the solution of the resulting ill-conditioned algebraic systems. An alternative approach recently proposed is the use of hyperbolic reformulations of dispersive models which allow for more efficient discretizations, [7,8,18,52].
Regarding high order methods, it is important to remark that while attaining high order in space is straightforward for DG methodologies, there are different possibilities concerning high order time discretizations. The original DG schemes of Cockburn and Shu employed high order Runge-Kutta schemes in time, leading to the family of RKDG schemes. An alternative consists in the family of fully implicit and semi-implicit space-time DG methods, see e.g. [19,82,98,99,[109][110][111]120,121]. Another different option that leads to high order explicit fully-discrete one-step schemes, and which is followed in this paper, combines ideas of the ADER approach of Toro and Titarev, originally developed within the finite volume framework [20,114,117,118], with space-time DG methods. This methodology, based on the ideas outlined in [41,43], makes use of an element-local space-time DG predictor, thus avoiding the cumbersome Cauchy-Kovalevskaya procedure of classical ADER schemes and thus allowing also the solution of complex PDEs in multiple space dimensions. Some examples of the wide range of applicability for this approach include the compressible Euler and Navier-Stokes equations, [40,41], compressible multi-phase flows [45], the Godunov-Peshkov-Romenski model of continuum mechanics, [17,47,93]. Discontinuous Galerkin schemes for hyperbolic PDE systems with non-conservative products have been proposed for the first time in [42,97], based on the ideas of path conservative schemes [22][23][24]85,89], which will be also a key point for the development of the numerical schemes proposed in this paper.
The rest of this paper is organized as follows: in Sect. 2 we present the original model [69] and a novel reformulation based on a decomposition of the specific Reynolds stress tensor P as P = QQ T . We furthermore introduce a viscous extension of the governing PDE system in order to define a rigorous and thermodynamically compatible vanishing viscosity limit of the model. We finally recall the Godunov formalism of thermodynamically compatible systems and prove that the proposed viscous system is thermodynamically compatible with the energy conservation law and with the entropy inequality. In Sect. 3 we present a novel thermodynamically compatible finite volume scheme, which mimics the aforementioned viscous extension of the system exactly at the semi-discrete level. In Sect. 4 a high order ADER discontinuous Galerkin method with a posteriori subcell limiter (MOOD) is presented for the new reformulation of the model proposed in this paper, including its viscous extension. Special care is taken concerning the conservation of total energy. Numerical results are shown in Sect. 5, where first a numerical convergence study is presented for third to sixth order ADER-DG schemes in space and time; subsequently, different schemes are compared with each other for three Riemann problems, discussing in particular the discretization of the nonconservative terms of the model in the context of thermodynamically compatible systems. The end of Sect. 5 contains numerical results for some challenging test problems for which experimental reference data are available, such as supercritical roll waves and the circular shock instability developing in the SWASI experiment, see [61]. The paper is rounded-off with some concluding remarks and an outlook to future work in Sect. 6.

Governing Equations
We consider the following overdetermined hyperbolic model for turbulent shear shallow water flows in multiple space dimensions, which has been recently proposed in [69] and which was also applied and studied in [11,25,80]: with the gravity constant g. The physical (primitive) state variables in (1)-(4) are the following: is the depth-averaged flow velocity and P = P(x, t) is the specific Reynolds stress tensor. For shallow water systems it is convenient to include the stationary bottom profile b(x) in the set of state variables. The reason is that this allows to represent stationary free surface waves associated with bottom jumps and to obtain well-balanced numerical schemes, see e.g. [21,22,64,89,90] for more details. Via straightforward calculations it can be shown that the system (1)-(4) admits the following extra conservation law with the total energy defined as h E = 1 2 gh 2 + hgb The bottom friction is taken into account by a coefficient C f and the dissipation function α is given according to [69] as

Reformulation of the Model in Terms of a New Evolution Variable
The above model requires trP ≥ 0 for hyperbolicity. In order to guarantee this property also at the discrete level for all times, we propose the following novel reformulation of the system (1)- (5). For this, we consider first the homogeneous part of equation (3) for the symmetric tensor P:Ṗ Here for shortness, for any f ,ḟ means the material time derivative:ḟ = f t + v · ∇, and L = ∂v ∂x = ∇v. Let us replace P by P = QQ T . What is the equation for Q? One obtains from (7): IfQ with an antisymmetric tensor B = −B T , the equation for P will be obviously satisfied. Thus, the equation for Q is defined up to an antisymmetric tensor B taking into account a proper rotation of the Reynolds tensor (for details, see [68]). We hypothesize that friction forces will drastically reduce the influence of this proper rotation, i.e. we take B = 0. Such a class of solutions is not equivalent to all solutions governed by the equation for P, but is able, as we will show, to describe complex flow configurations. What is a geometrical sense of such a decomposition P = QQ T ? Let us recall first the definition of the Gram matrix G (in the 2D case). Consider two vectors w i , i = 1, 2. The Gram matrix is defined as The 'dot' here is for the scalar product of vectors. It can be also written as i.e. the line vectors w i are the lines of Q. Let us recall that in our case P is the correlation tensor expressed in terms of the velocity pulsations (see [112] for details) as: Here the averaging operation, denoted by a "bar", is the depth averaging. The tensor P is positive definite due to the Cauchy-Schwarz inequality. Let us show that P can be presented in the form (10), i.e. there exist vectors w i : For this we take The last relation is well defined due to the Cauchy-Schwarz inequality: With P ik = Q im Q km written in terms of the new evolution variable Q and the notation ∂ m = ∂/∂ x m the above system can be rewritten again as an overdetermined PDE system as follows: with the conservative evolution variables h = h(x, t), hv = hv(x, t), Q = Q(x, t) and the stationary bottom It is easy to see that (3) is a consequence of (16) by simply multiplying (16) with Q T from the right and summing the transpose of (16) multiplied by Q from the left. It can be easily checked that also the new system (14)- (17) admits an extra energy conservation law which can be obtained as a consequence of (14)- (17). In terms of Q the total energy reads h E = 1 2 gh 2 + hgb + 1 2 hv i v i + 1 2 h Q i j Q i j , which for flat bottom b = 0 is a strictly convex function in the variables h, hv i and S i j = h Q i j . It is also a convex function of (h, hv i , Q i j ), if the turbulent energy is small compared to the gravitational potential energy (see "Appendix A" for details). Also note that due to trP = Q i j Q i j ≥ 0 the use of Q instead of P automatically guarantees a non-negative trace of P by construction, and hence also at the discrete level for all times. In this sense, system (14)-(18) is analogous to a so-called realizable turbulence model. At this point we emphasize that the thermodynamically compatible scheme proposed later in this paper will consider only the case of a flat bottom with b = 0.
Last but not least, we would like to point out the difference in the only apparently similar structure of PDE (16) and the governing PDE for the distortion field A ik in nonlinear hyperelasticity [47,93], which reads As one can easily see, the order of the matrix-product in the third term on the left hand side of (16) and (19) is exchanged. It is well-known that for hyperelasticity there is an additional conservation law associated with the determinant of the distortion field A ik and in the following we will show that the same applies to the determinant of the field Q ik . The time derivative of the determinant of Q can be easily obtained via the Jacobi formula, which expresses the derivatives of the determinant of a matrix in terms of the inverse of the matrix and the derivatives of the matrix itself: where Q −1 ki is a compact notation for (Q −1 ) ki . Applying (20) to (16) yields from which one obtains and therefore the sought additional balance law for the determinant |Q|, which for α = 0 has the same structure as the mass conservation equation (14). As such, we can assume that for h > 0 also |Q| > 0 holds. Via straightforward calculations it can be shown that for smooth solutions the conservation law (23) for the determinant |Q| is equivalent to the conservation law already found in [69]. Assuming α = 0 it reduces to

The Godunov Form of Nonlinear Systems of Hyperbolic Conservation Laws
In order to define the vanishing viscosity limit of system (14)- (17) and in order to introduce the new thermodynamically compatible finite volume schemes developed later in this paper, which are exactly compatible with the vanishing viscosity limit, it is necessary to recall the Godunov form [70] of hyperbolic PDE systems. We first consider only hyperbolic systems of conservation laws in two space dimensions of the type with flux tensor F = (f 1 , f 2 ), that admit the following parametrization according to Godunov [70] with the extra conservation law of the form where F k is the total energy flux in the k-th coordinate direction. Equations (29) and (30) are in the following called the Godunov form of the conservation law (28) and constitute an overdetermined system of PDE. The system is thermodynamically compatible if the following relations hold: Here, L is the so-called generating potential and E is the total energy density, which are the Legendre transforms of each other and thus satisfy We assume L and E to be strictly convex functions of their arguments, hence the transformation matrices between p and q variables, which are the Hessian matrices of L and E, respectively, verify It is easy to check that (30) is a consequence of (29), since scalar multiplication of (29) with p = E q yields which is the sought form of the total energy conservation law (30). For details on the class of symmetric hyperbolic and thermodynamically compatible (SHTC) systems and their application, see [17,47,70,71,[73][74][75]93,102]. The shallow water subsystem for flat bottom contained in (14)-(18) falls into the class of PDE (28)- (30). The corresponding potentials are and with the vectors q = (h, hv 1 , hv 2 ) T and p = (gh − 1 The associated Hessian matrices are and It is easy to see that with (40) and the flux tensor F = (hv k , hv i v k + 1 2 gh 2 δ ik ) T the energy fluxes (31) in (30) are which corresponds to the energy flux in (38).

Thermodynamically Compatible Vanishing Viscosity Limit
In order to define weak solutions for system (14)-(17), we define an associated thermodynamically compatible viscous system that satisfies at the same time an entropy-type inequality, as well as the total energy conservation law. In this section we assume a flat bottom with b = 0 for simplicity, as well as α = C f = 0, while a small parabolic dissipation term with dissipation coefficient ε > 0 is added to the equations. In order to guarantee exact total energy conservation, a non-negative production term T ik must be added to the governing PDE for Q: with the total energy Here, E 1 is the total energy potential of the shallow water subsystem (36)-(37) and E 2 is the total energy associated with the new object Q ik . In what follows, we will denote the inviscid part of the total energy flux in (47) by with the abbreviation that will be used later and which corresponds to the energy flux related to the shallow water subsystem, see also (43). The production term, T ik , which is needed to achieve the consistency of (44)-(46) with the total energy conservation law (47) reads The consistency with physics and experimental observations requires total energy conservation, see [68,69,80,100,101] for a more detailed discussion. In (50) the vector q = q i = (h, hv i , Q ik ) indicates the vector of primary state variables and E q i q j is the Hessian matrix of the total energy potential with respect to these state variables. One can show that the Hessian matrix is positive definite for small turbulent kinetic energy Q i j Q i j compared to gh, see "Appendix A" for details.

Theorem 1 (Energy conservation) The energy conservation law (47) is a consequence of equations (44)-(46).
Proof The shallow water subsystem (36)-(37) related to E 1 , which are the black terms in (44)- (45), directly falls into the general class of PDE (29)-(30) found by Godunov, hence the compatibility of the shallow water subsystem with the energy conservation law with energy potential E 1 is obvious. It is therefore enough to consider only the remaining terms associated with the quantity Q ik (red) and the viscous terms on the right hand side (blue). We first show compatibility of the red terms: (46) with the thermodynamic dual variables and considering only new contributions that are not yet contained in the Godunov-form yields After simple renaming of indices this proves the thermodynamic compatibility of the red terms contained in the left hand side of (44)- (46) with the red terms on the left hand side of the energy equation (47). We now consider the right hand side (blue terms): We define a viscous flux tensor g k as and a production term T that is equal to zero for all PDE apart from the non-zero production term T ik in the PDE for Q ik , see (46) and (50). Summation of the right hand sides of (44)- (46) with the thermodynamic dual variables p = E q yields The combination of the right hand sides of (44)-(46) therefore yields the right hand side of (47), which completes the proof.
Proof To see that the entropy inequality is a direct consequence of (55) we apply the Jacobi identity (20) to (55), which leads to from which one obtains With one obtains the following entropy-type inequality associated with system (44)- (46): where {i, j} ∈ {1, 2, 3}.
Throughout this paper, we will consider entropy solutions of (14)-(18) that satisfy (44)- (47) in the limit ε → 0. As shown later, the thermodynamically compatible scheme proposed in Sect. 3 of this paper is provably compatible with this vanishing viscosity limit, since it mimics the above viscous system exactly at the semi-discrete level. In the section containing the numerical results, we provide numerical evidence that also the high order ADER-DG schemes proposed in Sect. 4 of this paper as well as the numerical scheme already developed in [69] are compatible with this vanishing viscosity limit. The meaning of Theorem 2 is the following. The evolution equation for the tensor P (or for Q) is responsible for the vorticity transport, dissipation and production. While the transport and dissipation terms are clearly identified in previous works, it is not the same for the production terms. In the one-dimensional case the energy equation is equivalent to the 'entropy' equation and the 'entropy' (or vorticity) production is a consequence of the energy conservation. In the multi-dimensional case, the situation is completely different because the governing equations cannot be written in conservative form (the proof is given in [69]). So, the definition of weak solutions for such a non-conservative hyperbolic system which is compatible with the entropy production, should be given. In particular, such a definition was proposed in [69]. The Theorem 2 can be seen as a compatible alternative approach for the definition of weak solutions : the 'viscous' terms playing a major role in shocks guarantee the vorticity production. Moreover, the 'viscous' terms are consistent with the energy conservation law (Theorem 1) that is a necessary condition for all physically reasonable mathematical models.

Thermodynamically Compatible Finite Volume Scheme
In order to derive our new thermodynamically compatible finite volume scheme for system (14)-(18) that mimics the structure of the viscous system (44)-(47) exactly at the semidiscrete level, we proceed in a similar way as on the continuous level. First, a compatible scheme for the shallow water subsystem (36)-(38) is derived, based on a semi-discrete version of the Godunov form of (29)- (30). This corresponds to the discretization of the black terms in (44)- (47). Then, numerical viscosity together with an appropriate entropy production term is added to the scheme, which corresponds to the discrete analogue of the blue terms in (44)-(47). Last but not least the discretization of the red terms in (44)-(47) is discussed. To keep the presentation simple, we restrict ourselves to the one-dimensional case, but the generalization to multiple space dimensions is straightforward. To avoid confusion in the notation throughout this section we will use the lower case subscripts i, j, k, l, m for tensor indices and the lower case superscript r for the spatial discretization index. We emphasize again that the scheme proposed in this section is only valid for the flat bottom case with b = 0.

Compatible Schemes Without Dissipation Applied to the Godunov Form
A semi-discrete conservative finite volume scheme for system (28) in one space dimension based on the spatial control volume By adding and subtracting f r = f(q r ) we get We now try to obtain a discrete form of the total energy conservation law (30) also as a consequence of the discrete equations (62). For this purpose, we multiply (62) with p r = E q (q r ) from the left and get where the fluctuations D . We now compute the temporal rate of change of the sum of the total energy in cell r and r + 1, which yields It is clear that in order to obtain a flux conservative form of the discrete energy conservation equation we must require that the contribution of the left and the right fluctuation at the interface r + 1 2 is a flux difference, i.e.
where F r must be a consistent approximation of the total energy flux F. Inserting the definitions of the fluctuations into (65) yields Using the parametrization (29) and the associated relations (31) we get with F r = p r · f r − (vL) r and thus the sought relation that the numerical flux f r + 1 2 must satisfy is The condition (68) above is like a Roe-type property, but only for the vector f r + 1 2 instead of an entire Roe matrix. Based on the ideas of path-conservative schemes of Castro and Parés [22,89] we thus define the numerical flux via a path-integral in phase-space, since by the fundamental theorem of calculus we have for any path ψ = ψ(s) connecting p r with p r +1 , see also the pioneering work of Tadmor [108] for a similar construction of an entropy-conservative flux at the aid of a path integral. The last equality in (69) means a concrete parametrization of the chosen integration path using integration by substitution and a dimensionless integration parameter s in the range 0 ≤ s ≤ 1. In the following we choose two different parametrizations based on the simple straight-line segment path. Note that the choice of the path is arbitrary, hence we are free to choose a path that is somehow convenient for our purposes.
1. Segment path in the p variables (p-scheme). In the p-scheme, the path between p r and p r +1 is directly given by the straight line segment We thus obtain and therefore relation (69) results in By comparison with (68) we find that the thermodynamically compatible numerical flux of the p-scheme is therefore given by which by construction satisfies p r +1 − p r · f r + 1 2 p = (vL) r +1 − (vL) r and thus condition (68). The problem with the p-scheme is that it requires f in terms of p variables, which in general is very cumbersome, since usually f is easier known in terms of q rather than in terms of p. 2. Segment in the q variables (q-scheme). To avoid the above-mentioned problem, in the q-scheme the path between p r and p r +1 is now defined in terms of a straight line segment in the q variables, which means in terms of p variables the path is in general not a segment. We setψ Here, we only use the notationψ(s) to avoid confusion with the path used before in the p-scheme. We therefore have and thus condition (69) results in If we now check again condition (68) we still need to transform the jump in p variables into a jump in q variables. For that purpose, we define a Roe-type matrixL qq that satisfies the Roe propertyL which can be easily achieved by construction by the means of a path integral. In practice we first define the inverse of the Roe matrixL which is again a Roe matrix, but which is easy to compute, and which can be checked to satisfyẼ by construction. We thus obtaiñ which finally yields the desired thermodynamically compatible numerical flux of the q scheme as Note that if f is only easily known in terms of q variables, one can directly plug the straight segment path in terms of the q variables into the function f and into the Hessian E qq , without needing to compute p(q) at all! In practical calculations, we approximate all path integrals by numerical quadrature, which can be done up to any desired level of accuracy, see also [42,45,48,49] where this strategy has already been successfully used.

Compatible Scheme with Dissipation Applied to the Godunov Form
The above schemes are compatible with the parametrization (29) of the system (28) and also satisfy the extra conservation law (30). However, to obtain a dissipative scheme, we still need to add a compatible numerical dissipation. For that purpose we write a dissipative scheme for (28) of the form with the compatible flux f r + 1 2 as defined before and the additional dissipative numerical flux g r + 1 2 given by where μ r + 1 2 ≥ 0 is a scalar numerical dissipation. Henceforth we simply set with s r + 1 2 max an estimate for the maximum signal speed at the interface. For ϕ r + 1 2 = 0 this choice corresponds to a classical first order Rusanov-type scheme, see [107,115]. To reduce numerical dissipation in smooth regions, a TVD minbee flux limiter ϕ r + 1 2 is employed, which is defined as follows, see the second order TVD SLIC scheme described by Toro in [115], with the ratios of subsequent slopes of the total energy potential defined as Note that in regions of ϕ r + 1 2 = 1 the scheme exhibits no numerical viscosity at all. The production term T r will be defined later. Computing the dot product of (82) with p i yields where we denote the inviscid numerical flux for the total energy by F r + 1 2 . Since the dissipation-free scheme has already been shown to be compatible with the energy conservation law, in what follows, the explicit expression for F r + 1 2 is not needed, but is given here for completeness: We now rewrite the right hand side of (87) as The total energy flux including the dissipative terms thus reads as follows since the expression 1 2 (p r +1 + p r ) · Δq r + 1 2 is an approximation of the path integral using the simple trapezoidal quadrature rule. Making again use of the symmetric Roe matrix E r + 1 2 qq , which satisfiesẼ r + 1 2 qq (q r +1 − q r ) = p r +1 − p r , the semi-discrete total energy conservation law takes the form By requiring that we finally obtain the sought conservation form of the discrete total energy equation (87): The term T r is set identically to zero in all its components, apart from the equations that are needed for the entropy inequality, which are the nonconservative evolution equations for Q ik . Therefore, the term T r will be discussed later.
To summarize, the thermodynamically compatible dissipative numerical flux of the q scheme for (28) reads

Thermodynamically Compatible Discretization of the Terms Related to Q ik
We now present the discretization of the Reynolds stress tensor R ik = h P ik in the momentum equation that is thermodynamically compatible with the term (∂ m v i )Q mk in (46) and the term h P ik v k in the energy equation. To ease notation, we present the discretization only for the onedimensional case in x 1 direction. Extension to multiple space dimensions is straightforward. For the term (∂ m v i )Q mk we a priori choose the following discretization: Multiplication of the momentum equation (45)  : Using R ik = h Q im Q km and collecting terms leads to , we obtain the following compatible discretization for the numerical flux of the Reynolds stress tensor in x 1 direction: which needs to be added to the compatible dissipative flux (95) in the semi-discrete momentum equation, which then takes the form where f r + 1 2 hv i ,d is the part of the dissipative flux in x 1 direction (95) that refers to the momentum equation.
The last term in (46) that needs to be discretized is the convective term v m ∂ m Q ik , which requires compatibility with the mass conservation law (44) and the energy conservation (47). To achieve such a compatible discretization, the mass conservation equation needs to be multiplied with the remaining contribution E 2,h = E 2 and the PDE for Q ik is again multiplied with E Q ik , and the following condition must be satisfied: with the yet unknown average velocityṽ of the q-scheme according to the semi-discrete Godunov formalism, see (81). Collecting terms leads to from which we obtain the sought expression for the average velocity at the interface as v r + 1 In case the denominator is zero, we simply set the velocity to the arithmetic averageṽ r + 1 2 1 = 1 2 ṽ r 1 +ṽ r +1 1 . In order to get compatibility with the total energy conservation law also in the presence of numerical viscosity, we need to add the discrete production term to the PDEs of Q ik at the right and left element interface, according to the condition (93) already derived before: and The physical entropy production is always non-negative, since we assume μ r + 1 2 ≥ 0 and E r + 1 2 qq ≥ 0. It is obvious that (105) and (104) are discrete analogues of the continuous production term (50).
The final semi-discrete scheme for Q ik in one space dimension reads:

Summary of the Scheme and Stability Proof
For completeness, we now gather together all equations of the thermodynamically compatible scheme, thus obtaining with the fluctuations where f r q denotes the physical flux evaluated in cell r and f r + 1 2 q is the compatible flux for depth and momentum, i.e. for q ∈ {h, hv i }. Recall that the flux vector in the previous notation reads and is computed according to the q-scheme of the semi-discrete Godunov formalism presented previously. We have also introduced the fluctuations with R r + 1 2 i,1 being the compatible discretization of the Reynolds stress tensor in x 1 direction given in (99), and Let us also recall that the dissipative fluxes, g r ± 1 2 q , have been defined in (83).

Theorem 3 The semi-discrete scheme (107)-(109) admits the semi-discrete energy conservation equation
As a result the scheme is energy conserving and therefore marginally stable in the energy norm, i.e. the scheme satisfies Proof We first demonstrate that (114) is a direct consequence of (107)- (109). To this end we proceed as in the continuous case, i.e. we start by considering the time derivatives and we sum the contributions coming from equations (107) For the convective terms, the Reynolds stress tensor and the PDE related to Q we define Let us remark that the definitions of the fluctuations in (118) differ from the ones given in Sect. 3.1, where the terms related to the total energy E 2 associated with Q ik were not yet included. Finally, the dot product of p r by the vector of the blue terms in (107)-(109) yields Taking into account (104)-(105) in p r · T r = p r · T r − 1 2 ,+ + p r · T r + 1 2 ,− we get (93). Substitution of this result in (119) and making use of the developments in (89) gives Gathering (117), (120) and (118), we get (114): Let us now consider the discrete equation for the total energy, (114). Integrating it on a computational domain, Ω, we get Recalling that D r ± 1 2 ,∓ E represent the jumps of the energy flux at the interfaces and assuming the solution on the boundaries of Ω to tend to a constant value, the jumps of q are zero at the boundaries of Ω and then the boundary contributions of D E and also the contribution of the dissipation terms vanish. Hence, reordering the pairs of D r ± 1 2 ,∓ E to consider couples related to the interfaces instead of pairs corresponding to the cells yields Note that the summation over the blue dissipative fluxes is obviously a telescopic sum that vanishes. On the other hand, we can also prove that the contributions of the fluctuations at the interfaces reduce to a flux difference of the form To this end, we simply develop the fluctuations related to the total energy Reordering black terms and using (97) and (113) for the red ones, we obtain Finally, taking into account (101) with (103) in the above expression and the discretization of f r + 1 2 q given in (68), (81), we conclude which is the sought total energy flux difference. Therefore, under the hypothesis that the total energy fluxes on the boundary are zero, we have hence the scheme is marginally stable in the energy norm.

Path-Conservative ADER Discontinuous Galerkin Schemes
In this section we briefly recall ADER-DG schemes on rectangular equidistant Cartesian grids with a posteriori subcell finite volume limiter (SCL). The governing PDE system (14)- (17) can be cast into the following general form ∂q ∂t + ∇ · F(q, ∇q) + B(q) · ∇q = S(q, ∇q), the time, the state vector q ∈ Ω q ⊂ R m , the state space or phase-space Ω q ⊂ R m , the flux tensor F(q, ∇q) = (f, g), the nonconservative product B(q) · ∇q = B 1 (q)∂ x q + B 2 (q)∂ y q and the source term S(q, ∇q), which may also depend on gradients of the state vector. The general structure (129) is needed if we want to discretize also directly the underlying thermodynamically compatible viscous system (44)- (47), including the production term T ik . In order to solve very general nonlinear time-dependent PDE systems like (129) numerically, in this paper we employ the family of high order accurate fully-discrete path-conservative one-step ADER discontinuous Galerkin schemes supplemented with an a posteriori subcell finite volume limiter, see e.g. [17,41,42,50,56,124]. In the next sections we provide a brief description of the method and the reader is referred to the above references for more details. Concerning more details on the framework of a posteriori limiting (MOOD), the reader is referred to [29,38,39].

Unlimited High Order ADER-DG Schemes
The system (129) is discretized on a domain Ω making use of a uniform Cartesian grid with elements Here, x i = (x i , y i ) is the barycenter of Ω i and Δx and Δy are the the mesh spacings in the x and in the y direction, respectively. The numerical solution of (129) is defined in the space of piecewise polynomials of degree N and is denoted by u h (x, t n ). For each element Ω i it is sought under the form Here, ϕ l (x) = ϕ l 1 (ξ )ϕ l 2 (η) are the basis or ansatz functions, which are tensor products of one-dimensional ansatz functions ϕ l m (χ) on the unit reference element The mapping from the reference element to the physical one reads x = x i − 1 2 Δx + ξΔx and y = y i − 1 2 Δy + ηΔy with 0 ≤ ξ, η ≤ 1. The multi-index l = (l 1 , l 2 ) refers to the onedimensional basis functions ϕ l m that are employed in the tensor product. The basis functions on the reference element are defined as the Lagrange interpolation polynomials that pass through the Gauss-Legendre quadrature nodes of a Gaussian quadrature formula with N + 1 quadrature points. This choice automatically leads to an orthogonal nodal basis.
Multiplication of (129) by test functions ϕ k , which according to the Galerkin approach are chosen identical to the ansatz functions, and integration over Ω i × [t n , t n+1 ] leads to

Using (130) and integration by parts yields
where n is the outward-pointing unit normal vector at the cell boundary ∂Ω i , and q h is a local space-time predictor, the computation of which will be briefly explained later. Since in the DG framework the discrete solution is allowed to jump between two neighboring cells, a numerical flux is required on the boundary. For an exhaustive overview of numerical fluxes and Riemann solvers, see [115]. In this paper, we use the simple Rusanov-type flux with s max = max(|λ k (q − h )|, |λ k (q + h )|) + ε(2N + 1)/Δx being an estimate of the maximum signal speed at the interface, including also the viscous terms with viscosity coefficient ε, see [65]. In (133) q − h and q + h denote the boundary-extrapolated values of the space-time predictor from within the element and its neighbor, respectively. The non conservative products are discretized via a path conservative scheme, as forwarded by Castro, Parés and collaborators in [21][22][23][24]86,89] and which are based on the theory established in [85]. The term D q − h , q + h contains the jump in the non-conservative product and is computed at the aid of a path integral in phase space between the states q − h and q + h . Using the simple segment path the path integral reduces to The integral (135) is approximated via a simple trapezoidal quadrature rule. The use of path integrals based on the straight-line segment path is the common point between the path-conservative ADER-DG scheme presented here and the thermodynamically compatible semi-discrete finite volume method presented in the previous section.
Following [17,41,43,50] the predictor q h (x, t) is obtained at the aid of a weak formulation of (129) in space-time, which allows to completely avoid the Cauchy-Kovalewskaya procedure that was originally employed in ADER schemes, see [20,113,114,118,119].
The predictor solution is defined at the aid of space-time ansatz functions θ l = θ l (x, t) = ϕ l 0 (τ )ϕ l 1 (ξ )ϕ l 2 (η), which are again tensor products of the 1D basis functions ϕ l m (χ) and where now an additional temporal basis function is included, with t = t n + τ Δt: Multiplication of (129) by θ k and integration over Ω i × t n , t n+1 yields and after intergration by parts one obtains the final weak form in space-time: Equation (138) is a nonlinear element-local algebraic system in the unknownsq n l,i , while the coefficientsû n l,i are the known from u h (x, t n ) at the previous time. The solution of (138) is obtained by an iterative algorithm, whose convergence was proven in [17] for the case of hyperbolic conservation laws without non-conservative products and without source terms. Concerning the choice of a suitable initial guess for the unknown space-time coefficientsq n l,i , the reader is referred to [55,79]. This completes the description of the unlimited ADER-DG scheme.

A Posteriori Subcell Finite Volume Limiter
The numerical scheme presented above is high order accurate and linear in the sense of Godunov, hence it will inevitably generate spurious oscillations in the vicinity of shock waves and discontinuities according to the well-known Godunov theorem. In [12,46,50,124] a new a posteriori subcell limiter was introduced for ADER-DG schemes, using the ideas of the MOOD paradigm forwarded in [29,38,39] for finite volume schemes.
At the beginning of each time step, the unlimited scheme described in the previous section is run on the entire computational domain. This produces a so-called candidate solution, in the following denoted by u * h (x, t n+1 ). Next, the candidate solution is a posteriori checked against different numerical and physical detection criteria, such as the positivity of the water depth and of the determinant of Q. Furthermore, the absence of floating point errors (NaN) is required and we also require a discrete maximum principle (DMP) to be satisfied, see [50]. If any of these numerical or physical detection criteria is violated, a high order DG cell is marked as troubled and is scheduled for the a posteriori subcell finite volume limiting.
The cells Ω i that have been scheduled for subcell finite volume limiting are now split into (2N + 1) d finite volume subcells, which are denoted by Ω i,s with Ω i = s Ω i,s . This subdivision of a high order DG element into many small finite volume subcells does not reduce the time step of the DG scheme because the CFL number of explicit discontinuous Galerkin schemes scales with 1/ (2N + 1), while for the finite volume scheme used on the subgrid cells, the maximum Courant number allowed is of the order of unity. At time t n the numerical solution in the finite volume subcells Ω i,s is represented as usual via piecewise constant cell averages denoted byū n i,s and which are obtained from the high order DG polynomials u h (x, t n ) asū These subcell averages are now evolved in time at the aid of a second order MUSCL-Hancocktype TVD finite volume scheme with minmod limiter, which is also a predictor-corrector method and thus looks quite similar to the ADER-DG scheme. The main difference is that now the test function is unity, hence the volume integral over the flux term disappears, and the spatial control volumes Ω i are replaced by the sub-volumes Ω i,s : where the local space-time predictor q h is now easily obtained from the Cauchy-Kovalewskaya procedure, see [115]. Once the cell averagesū n+1 i,s of all subcells contained within cell Ω i have been computed at the new time t n+1 according to equation (140), the limited DG polynomial u h (x, t n+1 ) at time t n+1 can be simply obtained via a constrained least squares reconstruction. For this we require that and The constraint (142) means conservation of the solution within the element Ω i . In addition to the coefficientsû n+1 i,l of the limited DG polynomial, in all limited DG cells we also keep in memory the finite volume subcell averagesū n+1 i,s , since they serve as initial condition for the subcell finite volume limiter in the case when a cell is troubled also in the next time step, see [50]. This completes the description of the a posteriori subcell finite volume limiter. For more details, the reader is referred to [46,50,124].

Renormalization of Q
In order to maintain a strict compatibility of the discrete energy conservation law (18) with the discrete trace of P, for the inviscid case ε = 0 we proceed as follows: at the end of each time step, we compute in each degree of freedom of the DG scheme and in each control volume of the subcell finite volume limiter the trace of P from the total energy and subsequently rescale Q according to where the subscript l denotes a generic degree of freedom, Q n+1 l is the preliminary value as computed from the numerical scheme described previously andQ n+1 l is the final result of Q after rescaling at the end of each time step.
We stress that the above renormalization (143)-(144) is not performed for the case of a viscous system, i.e. when ε > 0, since for sufficiently fine meshes (well-resolved viscous flow), the compatibility with the energy conservation law (47) must hold automatically up to the order of accuracy of the numerical scheme, since Theorem 1 establishes the compatibility of (44)-(46) with (47) at the continuous level.

Numerical Tests
Throughout this section the gravity constant is set to g = 9.81. Moreover, we will use SI units : m, s, etc. without writing them explicitly.

Test Problem with Exact Solution
In the following we solve a test problem suggested in [69], which has an exact solution of the PDE system (1)-(5) that reads In our setup we choose the above exact solution at time t = 0 as initial condition for h and v, while Q is initialized as Q = diag( √ P 11 , √ P 22 ). Our numerical simulations are run until a final time of t = 1 in the computational domain Ω = [−1, 1] 2 with β = λ = γ = 0.1 and h 0 = 1. A third order ADER-DG scheme (N = 2) is run on a sequence of successively refined meshes of N x = 5, 10, 15, 20 elements. Since the exact solution is a global polynomial of degree one in space, the spatial discretization is exact for this test problem. Since time variations in this test problem are also quite small and since the high order DG schemes require a rather small time step for stability, the overall error that is observed on all meshes, see Fig. 1, is of the order of machine accuracy, as expected.
To complete the definition of the problem we set h 0 = 1, q 0 = 0.5 and C f = C r = 0. Let us remark that to get the sought solution, (147)-(148), a set of analytical source terms, calculated by substitution of (147)-(148) in (14)- (17), must be added to the right hand side of the original system. The simulation is run until t = 0.25 using ADER-DG schemes of polynomial degrees N ∈ {2, 3, 4, 5}. The errors in L 2 norm obtained for h, u and Q 11 are reported in Table 1. Overall, the expected order of accuracy is reached, see bold numbers in Table 1.

Riemann Problems
It is well-known that the numerical discretization of nonconservative hyperbolic PDE is notoriously difficult, see e.g. [4,23] for a more detailed discussion. This is particularly true when the nonconservative product is acting across genuinely nonlinear waves. This situation is usually the case in the nonconservative equation (16), which makes its numerical discretization particularly difficult. In [69] a special split scheme was developed, splitting the original system (1)-(5) into two quasi-conservative subsystems and during the solution of each of the subsystems, energy conservation was rigorously enforced. In this paper, two new unsplit schemes have been proposed for the discretization of (14)- (18). In the case of the path-conservative ADER-DG scheme described in Sect. 4, the total energy conservation law is explicitly discretized and for the inviscid case ε = 0 the object Q ik is renormalized at the end of each time step according to (143)-(144), in order to maintain discrete compatibility with total energy conservation for each degree of freedom of the DG scheme and for each subcell average in case of the subcell FV limiter. Instead, when the path-conservative ADER-DG scheme is applied to the viscous system with ε > 0 then no renormalization is carried out and the compatibility with the total energy conservation law must be guaranteed by the high order of accuracy of the scheme in combination with a sufficiently fine mesh alone, without any renormalization of Q ik . As such, the fully resolved direct numerical solution (DNS) of the viscous system (44)- (46), which also includes the production term, with small but not vanishing ε > 0, constitutes the highest level of fidelity concerning the discretization  of the non-conservative product since the solution can be considered as smooth and therefore there are no ambiguities concerning the proper definition of the non-conservative product at all. Instead, the compatible HTC scheme, which implements a semi-discrete Godunov formalism, does not directly discretize the energy equation at all, but total energy conservation is a mere consequence of all the other equations at the semi-discrete level. In this case, the discretization of the nonconservative products, of the Reynolds stress tensor and of the viscous terms (the red and blue terms in (44)- (46)) is carried out in such a manner that at the semi-discrete level total energy conservation is automatically ensured. In the following we solve three Riemann problems on the domain Ω = [0, 1] × [0, 0.5] with initial condition The left and right initial states, as well as the final simulation times t end , are summarized in Table 2. The values of the state variables not indicated in the table are set to zero, i.e. v 1 = 0, Q 12 = 0 and Q 21 = 0. The simulations are run with four different numerical schemes S1-S4: (S1) The split scheme of Gavrilyuk et al. [69], using a very fine mesh in one space dimension, which serves as a reference solution. This scheme makes explicit use of the total energy conservation law in each subsystem used in the splitting approach. (S2) A high order unsplit ADER-DG scheme (132) with polynomial approximation degree N = 3 and 11, 200 × 4 elements, applied to the viscous system (44)-(46) with small but positive viscosity parameter ε > 0 (vanishing viscosity limit). Since ε > 0 we do not solve the energy equation (47) explicitly and apply no renormalization to Q ik . To obtain the discrete compatibility with the energy conservation law (47) and for the proper definition of the nonconservative products, only a sufficiently fine mesh is needed in combination with the high order DG scheme (fully resolved DNS). This approach serves to generate an additional and totally independent reference solution. (S3) The new unsplit thermodynamically compatible HTC scheme based on the semidiscrete Godunov formalism described in Sect. 3. This scheme is by construction exactly compatible with the viscous system (44)-(46) at the semi-discrete level and therefore the semi-discrete energy conservation law is a direct consequence of the discretization of all the other equations and thus does not need to be discretized explicitly again. (S4) The high order unsplit ADER-DG scheme applied to the inviscid system, setting ε = 0 in (44)- (47). In this case, the energy equation (47) is explicitly discretized and the object Q ik is renormalized at the end of each timestep according to (143)-(144).
We emphasize that in all four schemes, discrete compatibility with the total energy conservation law is always assured in one way, or in another: either by directly using the total energy conservation equation within the numerical scheme (S1 and S4), or by achieving compatibility exactly at the discrete level (S3). In S2 the compatibility is merely achieved at the aid of negligible discretization errors by using a very high order scheme applied to the viscous system with ε > 0 and using a sufficiently fine mesh (fully resolved DNS). The computational results obtained for the three Riemann problems are shown in Figs. 2, 3 and 4, where also the exact mesh resolution is given for each scheme, together with the choice of the viscosity parameter ε in the case of the simulation of the vanishing viscosity limit. For all Riemann problems one can note an excellent agreement between the numerical solutions obtained with all four schemes (S1-S4) listed above.
This clearly highlights that discrete thermodynamic compatibility is a key feature for the correct discretization of nonconservative products that are acting across genuinely nonlinear fields, like the ones present in (16).

One Dimensional Brock Profile
Here we repeat the numerical experiment concerning one-dimensional roll waves proposed in [69] and compare the obtained numerical results with the experimental data provided by Brock in [15,16]. The initial condition is given according to [69] by h = h 0 (1+a sin(2π x/L) with a = 0.05, v 1 = gh 0 tan θ/C f , v 2 = 0 and Q = 1 2 ϕh 2 I. The bottom slope angle of this test problem is θ = 0.05011 with ∂ x b = tan(θ ), the still water depth is set to h 0 = 0.00798, the bottom friction coefficient is chosen as C f = 0.0036, the parameter C r is set to C r = 0.00035 and ϕ = 22.76, see [69]. The computational domain is Ω = [0, L] × [0, 0.5] with L = 1.3 and is discretized with 104 × 20 ADER-DG elements of polynomial approximation degree N = 3. Periodic boundary conditions are applied in x 1 and x 2 direction. Simulations are run for system (14)-(16) until a final time of t = 12.5. For this test, the bottom slope term is simply implemented as an algebraic source term in order to be compatible with the periodic boundary conditions. The numerical results obtained with the path-conservative ADER-DG scheme and the experimental profile of Brock are depicted in the left panel of Fig. 5. Overall, we can note a very good agreement between the numerical results and the experimental reference data. In the right panel of Fig. 5 a visualization of the a posteriori subcell limiter is shown (red cells are highlighted in red, while unlimited cells are plotted in blue). It can be noticed that the limiter is only active at the shock wave.

Numerical Simulation of the SWASI Experiment
In this last and most complex numerical test we carry out the simulation of the SWASI experiment proposed by Foglizzo et al. in [61] and which was numerically investigated in [80]. The flow field is turbulent and turbulence is responsible for the developing flow structures. The computational domain is Ω = [−1, +1] 2 and is discretized with a uniform Cartesian mesh composed of 256 × 256 ADER-DG elements of polynomial approximation degree N = 3.
The bottom topography of this test is given according to [80] by with r = x , L 1 = 0.02, A = 0.005, β = 0.07, R − = 0.08 and R 1 = √ 2. The model parameters for this test are set to C f = 0.0036, C r = 1 and ϕ = 2. The reference inflow discharge is chosen as q 0 = 1.2 · 10 −3 , while the reference water depth at R 1 is set to h 0 = 0.003. In contrast to [80], in this paper the initial velocity field is chosen to be the stationary equilibrium between bottom slope and bottom friction and which satisfies the following ODE in radial direction: with initial condition u r (R 1 ) = −q 0 /(h 0 R 1 ) for both regions, r ≤ R 1 and r > R 1 . This ODE is solved once at the beginning of the simulation at the aid of a classical fourth order Runge-Kutta scheme. Once the radial velocity u r is known, the water depth can be easily computed as h = q 0 /(ru r ) and the final velocity field is given by the equilibrium solution plus a sinusoidal perturbation as follows: with d = 0.005. We furthermore set Q = ϕh 2 0 I as initial condition for the object Q. The outflow through the central hole is generated by a sink term, setting h = 10 −2 , v = 0 and Q = 10 −5 I for r < 0.075 at all times.
The numerical simulations are carried out until a final time of t = 57. Figure 6 shows the computational results obtained for the water depth, the Froude number, the angular velocity and the instantaneous a posteriori subcell limiter map, in which red cells are highlighted in red, while unlimited cells are plotted in blue. The limiter is essentially activated along the moving shock front. The obtained numerical results agree qualitatively very well with experimental observations made in [61] and with the numerical results previously presented in [80]. In particular, one can note the characteristic cusp in the shock front that is visible in both, the experiment [61] as well as in the numerical results of [80]. At this point we would like to emphasize that the model formulation as well as the numerical scheme used in this paper are completely different compared to the model formulation and the scheme employed in [80]. While in [80] the problem was solved after rewriting the governing PDE system in polar coordinates, here we solve the problem directly in Cartesian coordinates on a Cartesian mesh. Furthermore, the model used in [80] was based on the temporal evolution equation of P, while here we use the new model formulation in terms of the object Q that guarantees trP ≥ 0 by construction. Last but not least, in [80] a split finite volume scheme was used, while the present paper employs an unsplit high order ADER-DG scheme. The fact that the numerical results obtained with different models and different schemes agree well with each other and with experimental observations shows the validity of the different mathematical model formulations as well as of the chosen numerical discretizations. As already pointed out in [80], the same simulation run with the same model parameters (C f = 0.0036) and the same numerical method on the same mesh applied to the simple shallow water equations leads only to a steady circular shock wave, without developing any shock instability and without showing the typical cusp of the SWASI experiment, see Fig. 7.

Conclusion
In this paper we have introduced a new reformulation of the first order hyperbolic model for unsteady turbulent shallow water flows introduced and studied in [11,69,80]. The main idea of the model reformulation proposed in this paper is the decomposition of the specific Reynolds stress tensor P at the aid of a new object Q so that P = QQ T . This guarantees that trP ≥ 0 by construction also at the discrete level for all times, since in terms of Q the trace of the Reynolds  7 Numerical simulation of the SWASI experiment with a fourth order ADER-DG scheme at time t = 57 s applied to the classical shallow water equations. Water depth (left) and Froude number (right). With the classical shallow water model no shock wave instability develops stress tensor, i.e. the turbulent kinetic energy, can be written as trP = Q i j Q i j ≥ 0. Compared to the previous model used in [11,69,80] we also add a thermodynamically compatible viscous flux and an associated entropy production term that together guarantee the compatibility of the viscous system with the total energy conservation law and with the entropy inequality, which in the new reformulation can be simply expressed in terms of an extra conservation law for the determinant of Q. Based on the Godunov form of hyperbolic conservation laws found by Godunov in his groundbreaking work An interesting class of quasilinear systems [70], we have derived a new thermodynamically compatible semi-discrete finite volume scheme that mimics the Godunov form of the inviscid conservative part of the system exactly at the semi-discrete level. The proposed schemes can therefore be called a discrete Godunov formalism, or a hyperbolic and thermodynamically compatible (HTC) finite volume scheme. Subsequently, also a thermodynamically compatible viscous extension of the scheme has been proposed, together with the thermodynamically compatible discretization of the remaining nonconservative terms and of the Reynolds stress tensor, which do not fit into the original Godunov formalism. At this point we stress again that the proposed scheme mimics the underlying viscous system of the mathematical model exactly at the semi-discrete level and as such also falls into the class of structure-preserving schemes, since all properties of the thermodynamic structure of the governing PDE system are properly maintained by the numerical scheme. The paper also considers high order path-conservative fully-discrete one-step ADER discontinuous Galerkin schemes with a posteriori subcell limiter that can be applied to both, the viscous and the inviscid form of the mathematical model. The performance and accuracy of all schemes is carefully assessed at the aid of three Riemann problems, where also a direct comparison with the scheme introduced in [69] has been shown. An excellent agreement between all different methods was observed in all cases. For the high order ADER-DG schemes a numerical convergence study was carried out at the aid of a manufactured solution, since the analytic solution used in [69] was too simple for a high order DG scheme. The new model was applied to the simulation of roll waves, comparing with the experimental data of Brock and obtaining an excellent level of agreement between the numerical and the experimental results. As a last test problem the we have carried out a numerical simulation of the SWASI experiment of Foglizzo et al. [61], using the computational setup proposed in [80]. Our simulations show the same cusp in the moving shock front that was already observed in the experiments and in the numerical simulations shown in [80].
In the future we plan to extend the new family of thermodynamically compatible schemes to the equations of nonlinear hyperelasticity [14,67,75,77,87,102] and to the unified hyperbolic model of continuum mechanics [13,17,47,93,102], as well as to hyperbolic reformulations of dispersive systems [7,18,37,58]. Further work will also concern the extension of the discrete Godunov formalism presented in this paper to higher order semi-discrete discontinuous Galerkin finite element schemes, see e.g. [36].
Another open challenge remains the development of thermodynamically compatible schemes like those presented in this paper that also maintain curl and divergence involution constraints exactly at the semi-discrete level, similar to the structure-preserving semi-implicit method recently proposed in [13], but which was not thermodynamically compatible.
an UniTN starting grant of the University of Trento. S.G. has been partially funded by the Excellence Initiative of Aix-Marseille UniversityA*Midex, a French Investissements d'Avenir programme AMX-19-IET-010, and the Russian Science Foundation (Project 20-11-20189).
Funding Open access funding provided by Università degli Studi di Trento within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Convexity of the Energy
Consider the total energy as a function of s = (h, r i = hv i , S i j = h Q i j ) T : We now use again Sylvester's criterion. The first three principal minors are positive. One can show that if the determinant of E q i q j is positive, the other principal minors are also positive.
It is equivalent to the inequality gh − Q 2 11 − Q 2 12 − Q 2 21 − Q 2 22 > 0. Thus the determinant is positive if the 'turbulent' energy is small compared to gh. In practice, it is always the case. In the following, for convenience, we consider the energy as a function of q = (h, hv i , Q i j ) T .