{'fraction_non_alphanumeric': 0.03296344647519582, 'fraction_numerical': 0.01747567054355566, 'mean_word_length': 4.1481594623491675, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A path conservative finite volume method for a shear shallow water model

The shear shallow water model provides a higher order approximation for shallow water flows by including the effect of vertical shear in the model. This model can be derived from the depth averaging process by including the second order velocity fluctuations which are neglected in the classical shallow water approximation. The resulting model has a non-conservative structure which resembles the 10-moment equations from gas dynamics. This structure facilitates the development of path conservative schemes and we construct HLL, 3-wave and 5-wave HLLC-type solvers. An explicit and semi-implicit MUSCL-Hancock type second order scheme is proposed for the time integration. Several test cases including roll waves show the performance of the proposed modeling and numerical strategy.


Introduction
The shallow water equations, sometimes also called the Saint-Venant equations, are used to model the flow of fluids in situations where the depth of fluid is small relative to the horizontal scale of the flow field variations [23]. They have been used to model the flow on the scale of the atmosphere and ocean, and have been applied for tsunami prediction, storm surges, flow around structures, etc. Because the model has only two independent spatial variables and does not require tracking the free surface, it provides a simpler approximation than the full three dimensional Euler equations with a free surface. The shallow water equations are derived from the incompressible Euler or Navier-Stokes equations by averaging them over the depth coordinate. The horizontal velocity is assumed to be weakly varying in the vertical coordinate which implies that the vertical shear is negligible. This allows us to ignore the second order velocity fluctuations and leads to a closed set of equations which we may be called the classical shallow water equations and are given by ∂h ∂t + ∇ · (hv) = 0 (1a) where h, v = (v 1 , v 2 ) are the water depth and velocity, b is the bottom topography, see Figure (1), g > 0 is the acceleration due to gravity, and we have included a frictional term to model the bottom friction with C f being the Chezy coefficient. The horizontal velocity v is a depth average of the three dimensional velocity field. Since the vertical shear is neglected, the classical shallow water equations cannot model large scale eddies ('roller') that appear near the surface and behind the hydraulic jump. Under the assumption of smallness of horizontal vorticity, a more general model called the shear shallow water (SSW) model can be derived [19,17,18,10] which includes the second order velocity fluctuation terms, and can be written as ∂h ∂t + ∇ · (hv) = 0 (2a) where the symmetric stress tensor P comes from the second order velocity fluctuations related to the perturbations from the depth average, and comprises of three independent components P 11 , P 12 , P 22 , and D is the dissipation tensor. Using Stokes type hypothesis, we can relate the dissipation tensor D linearly with the stress tensor P, and furthermore if we want the model to reduce to the classical shallow water model in the limit P = 0, we get [10] where the coefficient α is a function of the invariants of P. The equations (2) lead to an equation for "total energy" ∂ ∂t (he) + ∇ · hev + (gh 2 where e = 1 2 |v| 2 + 1 2 trace(P) + 1 2 Following [18,10], the coefficient α in the dissipation term is given by α = max 0, C r T − φh 2 T 2 , T = trace(P) = P 11 + P 22 Moreover, the quantities C f , C r , φ are model constants that must be calibrated using experiments. The equations for P would contain third order velocity fluctuations, see Appendix A, which can be ignored if the horizontal shear is weak [19] or modeled in some way so as to close the set of equations. The depth averaging process which is described in Appendix A reminds us of the well known Reynolds averaging of Navier-Stokes equations for turbulent flows which leads to a hierarchy of equations due to the non-linear nature and have to be closed with some turbulence model 1 . A similar model as above has been studied for the Favre-averaged compressible Navier-Stokes equations used to model turblent flows [3]. The SSW model (2) is non-conservative since the equations for P cannot be put in conservation form. The solution of non-conservative equations is a tricky issue when discontinuities arise since we need a proper notion of weak solution. The jump conditions will depend on the particular form of the equations and not all forms may yield the same jump conditions. The numerical solution of the set of equations (2) has been addressed in [10,4], both of which use a splitting approach. The approach in [10] splits the equations into two sub-systems called a-waves

Re-formulation of the SSW model
The SSW model (2) can be written in an almost conservative form. To do this, we define the symmetric tensors R ij := hP ij , The quantity R ij has units of stress while E ij has units of energy per unit volume. Then an elementary computation shows that the set of equations (2) can be written as the following set of non-conservative equations ∂U ∂t Note that the vectors B 1 (U ), B 2 (U ) are linear in U and in fact depend only on the momentum density m = hv It is usual to write the non-conservative terms in terms of a matrix-vector product, but since only the derivatives of h appear in the non-conservative terms of this model, it is more convenient to write it in terms of the vectors B 1 , B 2 as above. Coincidentally, this model is identical to the 10-moment model [15,2] from gas dynamics except for the presence of gravity terms and the non-conservative terms. While the set of equations (4) can be obtained by manipulating equations (2), the more fundamental way to obtain these equations is by the depth averaging process which is performed in Appendix A. It is also known that the quantity det P h is a convex entropy function [2]. Using the SSW equations, we can derive the entropy equation [10] ∂η ∂t The fact that the entropy is a convex function of U and this variable set arises naturally when we perform the depth averaging, indicates that the set of variables U has a more fundamental character than the set (h, m, P) used in (2). Hence we propose to construct numerical schemes starting from the SSW equations as given in (4). The numerical strategy we use is based on the concept of path conservative schemes [16] applied to system (4), and we recall the basic notions of path conservative schemes in a 1-D version of the above model in the next section. We will assume throughout the paper that the bottom topography b(x 1 , x 2 ) is a continuous function and independent of time, in which case the terms containing this quantity can be treated as source terms and need not be included in the Riemann solver.
3. Model in 1-D and notion of path conservative scheme We will consider the 1-D SSW model which can be written as where F = F 1 , B = B 1 and the source term is given by For simplicity of notation, we will sometimes write the velocity components as (u, v) = (v 1 , v 2 ). The system of equations (5) is a hyperbolic system with eigenvalues [10,3] The first and last eigenvalues correspond to genuinely non-linear characteristic fields in the sense of Lax [11], while the remaining eigenvalues correspond to linearly degenerate characteristic fields [10]. Hence λ 1 , λ 6 are associated with shock/rarefaction waves while the remaining eigenvalues give rise to shear/contact waves.
Ignoring the source term for the moment as they do not contain derivatives of U , let us write the non-conservative system as ∂U ∂t If we have discontinuous solutions, then we have to give a meaning to the derivative term which can be done by integration by parts if A is the gradient of a flux function as in case of conservation laws. If A is not the gradient of a flux, then the non-conservative product is interpreted as a Borel measure [8]. This definition requires the choice of a smooth path Ψ : [0, 1] × U × U → U connecting the two states U L , U R across the jump discontinuity at x = x 0 such that where U is the set of admissible states. Then the non-conservative product is defined as the Borel measure [8,12] µ( where δ is the Dirac delta function. This viewpoint is equivalent to the definition of non-conservative product proposed by Volpert [22]. Using this notion, a theory of weak solutions can be developed based on which the Riemann problem has usual structure as for conservative systems, leading to shocks or rarefaction waves corresponding to genuinely non-linear characteristic fields and contact waves corresponding to linearly degenerate fields. Across a point of discontinuity moving with speed S, a weak solution has to satisfy the generalized Rankinge-Hugoniot jump condition The choice of the correct path is a difficult question and has to be derived from a regularized model motivated from the physical background of the problem. In practice, it is usual to consider the linear path Then the jump condition for our model (5) becomes where The source term does not make any contribution to the jump conditions since it does not contain derivative of U . The Riemann problem is the building block of a finite volume method and this approach can be used for non-conservative systems also [12,16]. The main idea is to split the fluctuation into two parts corresponding to left moving and right moving waves arising in the Riemann solution, where the fluctuation is defined as The splitting of the fluctuation can be performed using a Roe-type Riemann solver or HLL-type Riemann solver, the latter being the approach taken in the present work. Assume that there are m waves in the Riemann solution with m − 1 intermediate states. Let us denote the wave speeds as S j , j = 1, . . . , m and the intermediate states as U * j , j = 1, . . . , m − 1 with U * 0 = U L and U * m = U R . The fluctuation splitting is given by Let us consider a partition of the domain into disjoint cells of size ∆x. Let U n j denote the approximation of the cell average value in the j'th cell at time t = t n . The first order scheme is given by For θ = 0 we obtain an explicit scheme and for θ = 1 we obtain a semi-implicit scheme; however the coupling in the semi-implicit scheme is only local to the cell. An exact solution process for the semiimplicit scheme is explained in the Appendix B. If the system is conservative, i.e., A = F (U ) for some F , then the above scheme can be written in conservation form with some numerical flux function [16].

Linear waves and jump conditions
In this section, we study the structure of the states forming the linearly degenerate waves. Let us define the average and jump operators by Then the jump conditions across a discontinuity moving with speed S lead to the following set of equations.
When we construct the HLLC solvers in later sections, we will use the information deduced in the following two sub-sections to decide the structure of the intermediate states.

Contact wave
Since the contact wave is linearly degenerate, the normal velocity u is continuous across the contact wave, u L = u R = u, and equal to the speed of the contact wave S = u. The jump conditions across the contact wave then yield the following four relations From the first and third conditions, we obtain Since we expect R 11 to be positive, the first factor cannot be zero and hence we require that v = 0, so that both velocity components are continuous across the middle wave. The second condition of (7) shows that R 12 is also continuous across the middle wave.
Remark. In the above derivation, we deduced that v is continuous across the contact wave but the jump conditions admit another solution in some special cases. If the two states are such that then v L = v R is an admissible set of states that satisfies all the jump conditions. The assumption R L 11 = R R 11 = 0 implies that there is no vertical shear which does not hold in practical situations that we are interested in. Hence it is reasonable to ignore this solution in the construction of the Riemann solver.

Shear waves
Let us now consider the simple waves corresponding to the eigenvalues λ 2 , λ 5 which are both linearly degenerate. If the two states U L , U R correspond to a λ 2 wave, then they lie on the same integral curve of the eigenvector corresponding to the eigenvalue λ 2 = u − √ P 11 . In terms of the variables (h, u, v, P 11 , P 12 , P 22 ), the corresponding eigenvector is [10] where c = √ P 11 , and the integral curve is given by We immediately see that h, u, P 11 are constant along the integral curve and such constants are also called Riemann invariants. This implies that λ 2 has the same value in the two states which is consistent with the fact that we have a linearly degenerate field. From the remaining equations, we can deduce that v √ P 11 + P 12 and det P are also invariant along the integral curve. Similarly, if we consider a λ 5 wave, the corresponding eigenvector is and the integral curve is given by We deduce that the quantities h, u, P 11 , v √ P 11 − P 12 , det P are invariant along the integral curve.
Remark. In contrast to the Euler equations, the contact and shear waves are here associated to different eigenvalues when R 11 > 0. In the limit when R 11 goes to zero, the contact and shear waves merge and the system loses its strict hyperbolicity since the eigenvectors r 2 , r 5 become parallel. Then the tangential velocity v is free to take any value along the integral curve. In order to capture the interactions between contact and shear wave, as well as the transition when R 11 becomes small, accurate numerical scheme should also considered in these waves. Indeed, interactions of contact and shear waves seems to play an important role for the considered model, as is observed in the two dimensional roll wave tests, where only the five wave solver which includes all waves in the Riemann solver is able to produce realistic looking solutions.

HLL Riemann solver
The HLL Riemann solver [13,9] includes only the slowest wave S L and fastest wave S R in the Riemann problem. There is an intermediate state U * between the two waves. The jump conditions (6) across the two waves are given by Eliminating F * we obtain the intermediate state The first three equations of (8) do not contain the non-conservative terms and we obtain the corresponding intermediate state The last three equations of (8) then yield Hence we have explicit solution for the full intermediate state. The split fluctuations are then obtained from We estimate the minimum and maximum speeds in the Riemann problem as follows which is similar to estimates used for Euler equations, see e.g., [9].

HLLC3 Riemann solver: 3 waves
The HLL solver does not include the linearly degenerate waves like the contact wave, which get excessively diffused in the numerical results. In the HLLC solver [20], the contact wave of speed S M is also included in the wave model so that we have three waves S L < S M < S R and two intermediate states U * L , U * R . We will determine the two intermediate states by satisfying the jump conditions across the three waves. Since the contact wave is linearly degenerate, the normal velocity in the two intermediate states is same and equal to the speed of the contact wave From the discussion on jump conditions in section (4.1) we know that v and R 12 are continuous across this wave so that v * L = v * R = v * and R * L 12 = R * R 12 = R * 12 . These states are shown in Table (1). We first consider the jump conditions across the S L , S R waves. The jump conditions for the continuity equation yield The jump conditions for the x-momentum equation yields Using (9), the above equation can be written as Substituting (10) into the first jump condition in (7), we can find an expression for the contact wave speed S M The jump conditions from the y momentum equation yield and since R 12 must be continuous across the contact wave, we obtain from R * L 12 = R * R 12 an expression for v * as We can observe that u * , v * are equal to the velocity of the intermediate state of the HLL solver. From the last three jump conditions corresponding to the energy equations, we obtain The intermediate states are finally given by and the fluctuations are given by

HLLC5 Riemann solver: 5-waves
The 3-wave HLLC solver improves upon the HLL solver by including the linearly degenerate contact wave. But the SSW model contains two more linearly degenerate waves and we can try to construct a multi-state HLL solver by including all linearly degenerate waves in the Riemann solution. When we include all the five waves, there are four intermediate states which are shown in Table (2). Note that the quantities h, u, R 11 are continuous across the two shear waves S * L , S * R as deduced in Section (4.2).
Jump across S L , S R waves. The jump condition across S L , S R waves is which can be written as The first jump condition gives h * α which is identical to equation (9). The second jump condition gives R * α 11 and it is identical to equation (10). The third and fifth conditions are coupled and their solution There is a common value of p in all the four intermediate states which is given by where the last equality can be derived from the h and x-momentum jump conditions across the S L , S R waves. Then the third and fifth jump conditions can be written as whose solution is given by From the above solution, we can compute R * α 12 = h * α P * α 12 and E * α 12 = 1 2 R * α 12 + 1 2 h * α u * v * α . The fourth jump condition gives E * α 11 which is identical to equation (13). The sixth jump condition gives E * α 22 The second and fourth wave speeds can now be estimated as Jump across u * wave. All the jump conditions are satisfied provided we choose u * as given in equation (11).
Jump across S * L wave. The non-conservative terms are absent for this wave since h is continuous. The first, second and fourth jump conditions are automatically satisfied. The third jump condition yields We can now check that the fifth jump condition is also satisfied.
Jump across S * R wave. The non-conservative terms are absent for this wave since h is continuous. The first, second and fourth jump conditions are automatically satisfied. The third jump condition yields We can now check that the fifth jump condition is also satisfied. All the intermediate quantities have been determined except v * * ; we have two estimates of R * * 12 from the jumps across S * L , S * R waves. Setting these equal to one another, we get an equation Finally, the fluctuations can be computed as 8. Second order scheme in 1-D We will follow a predictor-corrector approach like MUSCL-Hancock scheme which can also be thought of as an ADER scheme. We have to reconstruct the solution variables in each cell in order to get a better representation of the solution, and it is found to be beneficial in terms of getting non-oscillatory solutions to reconstruct primitive variables, which are taken to be Step 1. The first step is to predict the solution at half time level using a local formulation, i.e., without any information from neighbouring cells. We estimate the spatial variation of the solution by a limiter function, e.g., and transform to conserved variables The update equation is explicit if θ = 0 and is an implicit equation for U Step 2. The second step uses the predicted solution at half time level to update the solution to next time level

Two dimensional scheme
We consider a Cartesian mesh with cells of size (∆x, ∆y) where the cell centers are indexed as (j, k) and we use half indices (j + 1 2 , k) and (j, k + 1 2 ) to denote the cell faces. The MUSCL-Hancock type scheme can now be naturally extended to two dimensions by the following two step process.
Step 1. We estimate the the spatial variation of the primitive variables along the two directions and transform to conserved variables Now we can estimate the solution at the face centers by linear reconstruction, Then the solution in the cell is evolved by half a time step using the PDE (4) The above update is explicit if θ = 0 and is an implicit equation for U n+ 1 2 j,k if θ = 1 2 ; once we obtain U n+ 1 2 j,k , then we estimate the solution at the cell faces by a linear extrapolation from the cell centers Step 2. The second step uses the predicted solution at half time level to update the solution to next time level This completes the description of the higher order scheme in 2-D. The first order scheme in 2-D is easily obtained by eliminating the reconstruction step and the predictor step.

Numerical results
The methods developed in this paper are applied to some 1-D and 2-D test cases. The 1-D tests are performed using a purely 1-D code. In all the test cases, we take the acceleration due to gravity as g = 9.81 m/s 2 . The time step is chosen from a CFL condition of the form ∆t = CFL 1 max j,k λx(U j,k ) ∆x + λy(U j,k ) ∆y where λ x , λ y are the maximum wave speeds along the x, y directrions respectively, and we use CFL = 0.5 in all test cases. During each computation, we monitor positivity of the values of h, P 11 , P 22 and we do not perform any artificial clipping of negative values.

1-D shear test problem
This is a Riemann problem without source terms which gives rise to two shear waves. The initial depth, normal velocity and the stress tensor are constant in space and are given by  [4] on a grid of 10000 cells. The solution consists of two shear waves with discontinuities only in the transverse velocity, the P 12 and the P 22 components of the stress tensor. It is not surprising that the numerical error is larger for solvers that do not have these waves in their intrinsic structure (HLL and HLLC3). However, the differences are reduced when second order accurate schemes are used and globally the mesh convergence is achieved for all schemes and are in accordance with the solution obtained in [10], [4].

1-D dam break problem
The dam break problem models a situation where a dam gate is suddenly opened. The initial condition has a jump only for the depth variable and all the remaining variables are constant in space. We solve the Riemann problem without source terms in which the initial velocity is zero everywhere and the components of the stress tensor are P 11 = P 22 = 10 −4 m 2 /s 2 , P 12 = 0 m 2 /s 2 , while the initial depth has a discontinuity given by The solution contains one rarefaction wave and one shock wave separated by a contact discontinuity. In this context, the stress tensor is initially very small and almost at the limit where there is no shear and det(P) is also small (here it is 10 −8 ). The first order results are shown in Figure (4) and second order results are shown in Figure (5). The numerical solutions show some differences depending on the Riemann solver used. The HLLC3 and the HLLC5 solvers converge asymptotically to the same solution. However, the HLL solver converges asymptotically to a slightly different solution. This is clearly visible in the shape of P 11 . Moreover, none of these solutions matches with the solution obtained in previous works [4]. Such differences can be expected since we are using different jump conditions for the non-conservative system. There are also differences between the Riemann solvers developed in this work. The HLLC3 and HLLC5 solvers converge to the same solution under grid refinement while the HLL solver shows some differences compared to these two. In this sense, the inclusion of the intermediate linear waves in the Riemann solver seems to be important.

1-D modified dam break problem
This is a modified dam break problem where additional jump is initially introduced in the transverse velocity component. The Riemann problem is solved without any source terms and the initial condition is given by while the remaining quantities are constant in space and are given by v 1 = 0.1, P 11 = P 22 = 4 × 10 −2 , P 12 = 10 −8 . The solution shown in Figure (6) is obtained from the second order scheme and behaves similar to the dam break problem with two additional contact discontinuities. However, contact discontinuities move slowly than for the dam break problem and det(P) is less close to zero (of order 10 −3 ). With this modified problem, all the Riemann solvers asymptotically converge to the same solution. As expected, the resolution is improved when the approximate solver contains more physical waves.

1-D roll wave problem
This problem models the flow of a thin layer of liquid down an inclined bottom and we include all the source terms. The initial condition is given by Here θ = 0.05011 rad is the inclination angle, C f = 0.0036, h 0 = 7.98 × 10 −3 m, a = 0.05, φ = 22.76s −2 , C r = 0.00035, L x = 1.3 m. The bottom topography is given by b = −x tan θ and the boundary conditions are periodic. The chosen parameters correspond to those considered in [14] which leads to the formation of 1-D roll waves starting from a uniform flow which has the same structure as in Brock's experiments [6,7]. The water depth is compared with Brock's experimental data in Figure (7), which shows a hydraulic jump around x/L x = 1 but also a smooth profile immediately behind it which is the signature of a roll wave. The classical shallow water model does not predict this roll wave profile but only gives rise to the hydraulic jump. Figure

2-D analytical solution
The rate of convergence can be estimated from a 2-D exact solution developed in [10] which is linear in space and non-linear in time. The exact solution is given by where the parameters are taken to be h 0 = 1 m, λ = 0.1 m 2 /s 2 , γ = 0.01 m 2 /s 2 , β = 10 −3 /s. Previous works [10], [4] have shown first order convergence of the error norm for this test case. We perform the computations on the domain [0, 10] × [0, 10] upto the time t = 50 seconds using Dirichlet boundary conditions where the solution in ghost cells is set to the exact solution. The HLL solver is run on meshes 10 2 , 20 2 , 40 2 , 80 2 , 160 2 while HLLC3 and HLLC5 solvers are run on meshes 10 2 , 20 2 , 40 2 , 80 2 . The convergence of the error norm is shown in Figures (9) for all three Riemann solvers, which shows second order convergence of our method. The HLL solver requires finer meshes before we see the second order convergence probably due to its larger numerical dissipation as a consequence of not resolving the linear waves.

2-D roll wave problem
This problem is a 2-D extension of the 1-D roll wave problem described in Section (10.4) and models the flow of a thin layer of liquid down an inclined plane. The initial conditions are given by The chosen parameters correspond to those considered in [14] which leads to the formation of 1-D roll waves starting from a uniform flow which has the same structure as in Brock's experiments [6,7]. A small two dimensional perturbation is added to the water depth whose amplitude is controlled by the coefficient a in the initial condition. The boundary conditions are periodic in both directions. Figure (10) shows the surface plots of the depth field on a 1040 × 400 mesh while Figure (11) shows similar plots on a finer mesh of 2080 × 800 cells. The first order scheme shows rather smooth solutions due to higher numerical dissipation while second order scheme shows more features in the solution. The HLL and HLLC3 solvers also show somewhat smooth solutions while HLLC5 shows a more turbulent-type of solution. However, even the HLL/HLLC3 solvers show transverse wave structures in the second order scheme. The first order results from HLLC5 solver look similar to the corresponding results in the literature [10], [4]. The development of the water depth profile with time using HLLC5 solver can be seen in Figure (12), which shows three distinct phases. In the first phase, the profile develops a 1-D structure similar to the 1-D computation in Section (10.4) and is similar to Brock's experimental profile. At around time t = 10, the profile develops some perturbations at the hydraulic jump which generates transverse structures and subsequently this spreads to the whole domain. In Figure ( Mesh size ∆x 10 −10 10 −9 Mesh size ∆x 10 −10 10 −9 (c) HLLC5  ; the HLLC3 solver yields somewhat regular wave patterns while the HLLC5 solver shows a more dis-organized behaviour in the solution. We compute the spectrum of kinetic energy of the velocity fluctuations which is shown in Figure (15). Firstly, the spectra at times t = 36 and t = 60 are very similar which indicates that a statistically steady state has been reached. All the three solvers show a k −4 spectrum in intermediate wave-numbers. This is in contrast to two dimensional hydrodynamic turbulence which exhibits a k −3 energy spectrum [5]. We also show the k −5/3 line which is observed in three dimensional turbulence. With the present computations, we are not able to definitely conclude on the structure of the energy spectrum in the small and intermediate wave-numbers and a more refined computation with a higher order scheme may be necessary to clearly identify the precise scaling law in those wave-numbers. These results show that a higher order scheme is essential to capture the turbulent-like structure in the solution. Also, there seems to be a critical dependence on the wave structure included in the Riemann solver and the type of solutions observed in the computations. Only the five wave solver that models all the waves and combined with a higher order scheme yields solutions which seem to correspond better with observations in experiments [1].

Summary and conclusions
The present work deals with the shear shallow water model which is a higher order version of the classical shallow water model since it includes effects of vertical shear. The model has a non-conservative structure which poses difficulties in its numerical solution. Previous works have developed Riemann solvers by splitting the model into two parts, while in this work, we develop a unified Riemann solver using a more fundamental form of the equations which arises directly from the depth averaging process and is similar to the 10-moment equations in gas dynamics. By using a linear path in conservative variables, we develop path conservative Riemann solvers for this model. In particular, we develop HLL, 3-wave and 5-wave HLLC solvers and apply them to a set of test problems in 1-D and 2-D. A second order version of the scheme is developed using a predictor-corrector approach and the order of accuracy is checked numerically, while existing works have only demonstrated first order results. Among the three Riemann solvers, the results demonstrate that the HLLC5 solver is the most accurate in terms of resolving all the waves that arise in the solution. The 1-D results compare well with results from previous 5-wave solver available in the literature except in the 1-D dam break problem, where we see some differences in the shock. Such differences can be expected due to different jump conditions in our model compared to previous approaches. The 1-D roll wave problem develops the roll wave profile and compares well with Brock's experimental results. The 2-D roll wave problem shows significant differences in the performance of the Riemann solvers. Only the HLLC5 solver develops a turbulent-like solution which is similar to previous results in the literature and such a solution seems to correspond to real flows, see e.g. [1]. The solution has a well defined y-average with fluctuations super-imposed on top of this. Moreover, the fluctuations have a kinetic energy spectrum with well-defined scaling laws being visible, though more refined calculations are necessary to make definite statements on the precise structure of the solution. The higher order scheme using the 5-wave solver developed in this work is thus very promising for further studies on shear shallow water flows.

Appendix A. Derivation of the SSW model
The gravitational force is assumed to act in the negative x 3 direction. Let x 3 = ξ(x 1 , x 2 , t) denote the height of the free-surface and x 3 = b(x 1 , x 2 , t) denote the location of the bottom surface, see Figure (1). Then the depth of water is h = ξ − b. The flow is assumed to be governed by the incompressible Euler equations which are given by 2 where (u 1 , u 2 , u 3 ) is the 3-D velocity field, p is the pressure and we assume that the density ρ is constant. A material point on the free surface F (x 1 , x 2 , x 3 , t) = x 3 − ξ(x 1 , ξ 2 , t) = 0 remains on the free surface which means that 3 dF dt = 0 =⇒ ∂ξ ∂t + u α ∂ξ ∂x α − u 3 = 0 on  Define the depth average φ = φ(x 1 , x 2 , t) of any quantity φ = φ(x 1 , x 2 , x 3 , t) by The fluctuation with respect to the average value is φ = φ − φ and clearly φ = 0 Moreover, we note the following identities Integrating the continuity equation over the depth of water and using (A.3) yields where D/Dt denotes material derivative, the prime quantities are non-dimensional and Fr = U/ √ gH is Froude number. Ignoring terms of O(ε 2 ), we obtain the hydrostatic approximation where p a is the atmospheric pressure at the free surface which may be taken to be constant. Hence IfT ≤ φh 2 then the solution isT itself; otherwise we solve the quadratic equation f (T ) = 0 and choose the unique positive root in (φh 2 , ∞) as the solution. Once T is obtained, we get P 11 , P 22 from (B.1), (B.2); and also P 12 can be determined from the E 12 equation which is of the form 1 2 hP 12 + α|v| 3 θ∆tP 12 = S 12 :=Ũ n+1 Since we know the tensor P, we can update the energy tensor E which completes the update of all the quantities.