Shallow Moments to Capture Vertical Structure in Open Curved Shallow Flow

Abstract Shallow water models are successfully used for simulating geophysical flows like river floods, tsunamis, sediment transport, or debris flows. Depth-averaged models are in general attractive due to their low computational cost. However, information on the vertical velocity is required a priori, typically by assuming a parametrized profile using a time-independent function, e.g., a constant. A systematic generalization of depth-averaged flow models is given by the shallow moment method. This method retains transient information about the vertical flow profile by using a finite Legendre expansion for said profile with time and space-dependent coefficients. The shallow moment approach allows to include vertical information and generates a hierarchy of equations that in the limit, recover the reference system before depth-averaging. Even a low number of basis functions significantly increases the model’s capability to resolve vertical information and, whenever this is relevant for the modeling task at stake, significantly improves the predictive power compared to classical shallow water systems. In the paper, we show a two-dimensional simulation result based on the shallow moment system to predict processes unobservable in classical shallow flow, namely secondary flow and the redistribution of velocity profiles, which are of relevance in curved channelized flow as experimentally demonstrated by Steffler.


Introduction
Depth-averaged flow equations are a central element of geoscientific mathematical models describing global-, regional, and local-scale geophysical phenomena like atmospheric flow (Kalnay 2002), tsunamis (Gisler 2020), river hydraulics and channel flow (Yotsukura and Sayre 1976), as well as geomorphological surface flow such as snow avalanches (Hutter, Wang, and Pudasaini 2005), rockfall or debris flow type landslides (Kowalski and McElwaine 2013).The first and also most prominent example of a depthaveraged mathematical free-surface flow model is the classical shallow water system, also referred to as the St. Venant equations in one space dimension.For most realistic situations, however, the underlying assumptions of the shallow flow theory are too idealized.As a consequence, the scientific community has derived an overwhelmingly large number of adaptations to the original theory to account for flow at high internal friction, strong coupling with the bottom topography, or multi-component and multi-phase behavior.In contrast to the classical shallow water equations, which assume a constant vertical velocity profile, some depth-averaged shallow flow systems establish quasi-static, yet material-dependent or heuristically motivated vertical velocity profiles by a so-called shape factor (Pudasaini and Hutter 2007).Separating non-hydrostatic pressure effects furthermore allows to capture dispersion -a well-known mathematical model for this type being the Green and Naghdi (1976) model.In a more recent approach, the shallow moment method (Kowalski and Torrilhon 2019) has been introduced to account for transient changes in the velocity profile based on a Legendre expansion of said profile.
A regularized version of the shallow moment method was proposed in Koellermeier and Rominger (2020), circumventing the potential loss of hyperbolicity for higher-order systems.The latter has been proven useful to simulate geophysical processes, such as sediment transport (Garres-D� ıaz et al. 2020).The physical regime considered in the present work, however, did not suffer from a loss of hyperbolicity.We therefore performed our analysis on the non-regularized model.
Retaining information about the vertical flow structure, which is otherwise lost in classical depth-averaged shallow flow models, is for instance relevant in physical shallow flow regimes where frictional losses and mass in-take at the bottom topography depend on the flow's strain rate, hence the flow's vertical velocity profile.Significant strain rates naturally arise for instance in stratified density flows (Kowalski and McElwaine 2013) or sediment transport (Garres-D� ıaz et al. 2020).An accurate representation of strain rate information at the bottom of the flow requires a significant number of moments to be considered.Even a low-order shallow moment system, however, provides information on the gradient that allows the construction of new closure relations.This idea is not new and has for example also been investigated in Zhang et al. (2013) for the Boussinesq-Green-Naghdi equations.In contrast to the latter, however, the here proposed shallow moment approach provides a strictly hierarchical system structure that can be systematically extended to higher modes of the resolved velocity profile.This allows for the construction of numerical solvers that can traverse the moment hierarchy, hence enabling the study of different models within the same numerical framework.In contrast to this, comparing different flow models typically requires new, tailored numerical solutions on a case-by-case basis.
The relevance of vertical information even within shallow flow settings is for instance shown in curved channel flow experiments, such as one presented by Steffler (1984).The experiment considers a circular, channelized shallow flow situation and is particularly suited for our purpose as it provides measurements on the flow profile in the vertical direction.This provides insight into four interesting flow effects, see Figure 1: � Secondary flow, namely a circulation in cross-flow direction.� Superelevation, a water surface gradient implying a higher water level at the outer channel boundary than the inner boundary due to centrifugal forces.� Redistributing velocity profiles, denoting a changing vertical velocity profile at different positions of the channel due to different flow conditions.� Transient velocity profiles, signifying a temporary evolving vertical velocity due to different flow conditions.
While superelevation effects can at least be qualitatively captured via classical shallow flow theory, secondary flow as well as transient and redistribution effects cannot be represented at all.Rather different compensation strategies exist that introduce shape factors for correction.
For the two-dimensional shallow moment system, numerical simulations are currently only available for the axisymmetric system (Verbiest and Koellermeier 2023), using a relaxed formulation of the shallow moment approach proposed in Koellermeier and Rominger (2020).In this article, we present two-dimensional results in a Cartesian reference system.The mathematical model has the structure of a (conditionally) hyperbolic system with non-conservative contributions and source terms.We outline the structure of our numerical solver based on the theory provided in Castro (2006Castro ( , 2009aCastro ( , 2009b;;Canestrelli et al. 2009), as well as the references contained therein.While our implementation is currently restricted to first order, the mathematical model is suited to be solved within an ADER-WENO-type framework as proposed in the references mentioned above.
The paper is organized as follows: We start by briefly introducing the reference system used for the shallow moment approach in section 2 based on a simplified version of the incompressible Navier-Stokes equations.In section 3, we introduce the shallow moment approach in two dimensions, followed by section 4 where we compare our model to the model derived by Ghamry and Steffler, called Vertically Averaged Momentum (VAM) model (Ghamry and Steffler 2002).The VAM equations constitute a dispersive model also based on expanding the horizontal velocities as well as the pressure.
In section 5, we discuss suitable finite volume schemes to solve the governing equations.In section 6, we present Steffler's experiment, a laboratory experiment studying the flow dynamics of curved channels, and present our numerical results.We conclude with a summary and outlook in section 7.

Reference system
The governing equations are given by the incompressible Navier-Stokes system, where the flow is bounded by a basal topography h b ðt, xÞ and the free surface h s ðt, xÞ: Following (Kowalski and Torrilhon 2019), we assume kinematic boundary conditions (Pudasaini and Hutter 2007), conduct a scaling analysis based on a shallow flow approximation, and assume a hydrostatic pressure-dominated flow, before introducing a mapping of the z-axis onto the unit interval f 2 ½0, 1� using the transformation As a result, we obtain the still vertically resolved reference system, in which ~indicates mapped quantities Here, h(t, x) denotes the height of the flow, namely the difference between free surface and bottom topography, ðũðt, x, fÞ, ṽðt, x, fÞÞ stands for the velocity in x-and y-direction respectively, and ðr xz ðt, x, zÞ, ryz ðt, x, zÞÞ are the components of the deviatoric stress tensor.The density q and the gravitational acceleration g i ¼ ðge x , ge y , ge z Þ are constants of the system.
We furthermore introduced the mean velocities u m ðt, xÞ :¼ and the vertical coupling operator Notice that all vertical flow information is concentrated in the coupling operator.In the absence of vertical coupling, e.g., when the vertical velocity profile is constant ðũðt, x, fÞ ¼ u m ðt, xÞÞ and ðṽðt, x, fÞ ¼ v m ðt, xÞÞ, the classical shallow water model can be recovered by choosing an appropriate material closure relation for r xz and r yz .
Additionally, we augment the existing reference system with the trivial equation for the bottom topography for numerical reasons discussed later in section 5.However, the equation can be extended to model a transient bottom topography by including appropriate fluxes.The Shallow-Moment-Exner model introduced by Garres-D� ıaz et al. (2020) demonstrates this in a one-dimensional setting.
The momentum balance in z-direction simplifies to the hydrostatic pressure relation p ¼ hð1 − fÞqge z and is already substituted in (2) and (3).Therefore, no direct evolution equation for the velocity w in z-direction is available.However, by integration of the mass balance in its form as a divergence constraint hr � ðu, v, wÞ T ¼ 0 multiplied by h, we can recover an explicit expression for the vertical velocity w :

Bottom and free surface friction
To close for the stress tensor, we impose boundary conditions at the surface h s and the bottom topography h b of the flow.In many geophysical applications, including the curved channel, the flow is highly influenced by the bottom friction and it is sufficient to treat the surface h s as a stress-free boundary For the bottom friction, the original paper (Kowalski and Torrilhon 2019) states a slip boundary condition given by where k 0 ½m� is the slip-length and l 0 N s m 2 � � is a characteristic viscosity, introduced such that the slip length k 0 has unit ½m�: Other choices for the boundary conditions are also possible.Alternatives were for instance introduced by the authors of Garres-D� ıaz et al. (2020), who used the Manning friction law to describe a sediment transport model based on the shallow moment approach.Very similar to their approach, we use the Chezy friction model as used in Vasquez, Millar, and Steffler (2011) and is stated as The dimensionless Chezy coefficient C � is defined by a relation depending on the fluid height h and the effective roughness height k s :

Bulk friction
The simple and well-suited bulk friction model for water flow is described by a Newtonian fluid with constant viscosity l 2 R þ : The closure for the deviatoric part of the stress tensor is then given by Note that the factor h −1 is a consequence of mapping the z-axis onto the unit interval.

Shallow moment system
In this section, we derive the shallow moment system using an orthogonal Legendre basis (on the interval I ¼ ½0, 1�).We start by deriving the general result of projecting (1)-( 3) followed by some discussion on the projection of the material closure.Lastly, we rewrite the system in its abstract matrix form.A derivation for more general types of basis functions is possible.A onedimensional derivation for Chebyshev polynomials and non-orthogonal basis functions can be found in Steldermann, Torrilhon, and Kowalski (2022).
The fundamental idea of the shallow moment method is to approximate the velocities ũ and ṽ by using a finite expansion of order N (from now on called level and the components are indicated by the index k) with an appropriate set of basis functions / i ðfÞ and coefficients a i ðt, xÞ In the above notation, we used the fact that / 0 ¼ 1 for Legendre polynomials and therefore The mean velocity is an important quantity in order to recover the classical shallow water equations on the zero level N ¼ 0 and we may use the expressions u m and a 0 (or v m and b 0 ) interchangeably, depending on the context.In principle, we can use different levels for a i and b i , however, since we do not assume a preferred flow direction, both levels are considered equal.
The fully resolved velocity fields ũ and ṽ are now replaced by its ansatz using ( 13) and ( 14).The unknown coefficients a i and b i can be obtained by using a Galerkin projection of the momentum balance in both, the x and y direction with respect to each basis function / k h�, / k i :¼ Computing the moments of the momentum balances (2) and ( 3) without yet inserting the closure relations for the stress tensor yields the system ha i a j h/ i / j , / k i |ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl } ha i b h/ i / j , / k i |ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl } and ha i b j h/ i / j , / k i |ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl } Note that the first two terms of the right-hand side of ( 15) and ( 16) correspond to the Galerkin projection of the vertical coupling operator x and used the integral of the basis function / k defined as U k ðfÞ :¼ Ð f 0 / k ðfÞdf: Partial integration of the stress tensor leads to the following two terms on the right-hand side, where the model-dependent closure relations for the bulk stress and friction at the boundaries can be inserted.
Since we consider scaled Legendre polynomials orthogonal on the unit interval, the mass matrix M ik is also orthogonal.Therefore, inverting the mass matrix reduces to a simple division by its diagonal elements.Additionally, the vector W k ¼ ð1, 0, :::, 0Þ T , since / 0 ¼ 1 is part of the basis and therefore orthogonal to all other basis functions.
In matrix notation, the general structure of the shallow moment equations ( 1), ( 15) and ( 16) can be written as a set of evolution equations with non-conservative products and a source term with with k ¼ 0, :::, N and an already used the stress-free, free-surface boundary condition (7).

Projection of the bulk friction
For a Newtonian fluid as described in ( 11), the Galerkin projection yields and Note that the bulk friction term is only present for N � 1: For the shallow water or equivalently the shallow moment model with N ¼ 0, friction can therefore only be prescribed with a suitable bottom friction closure.

Projection of the bottom friction
To ease the notation, it is convenient to use a normalized Legendre basis with / k ð0Þ ¼ 1 8k: Choosing normalized basis functions do not influence the results derived above.
For the slip friction model of (8), we obtain and Projecting the Chezy friction model (9) yields ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi and ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi In conclusion, while a Newtonian constitutive relation is certainly the most common closure for a fully resolved, three-dimensional incompressible flow model, it cannot be accounted for as a bulk friction term in the classical depth-averaged shallow water limit (N ¼ 0).The shallow moment hierarchy, however, allows to inherit closure relations from fully resolved models.In particular, they rely on the availability of gradients of the vertical velocities, allowing closure relations of the form.rðh, u, r � u; hÞ, where h are the parameters of the closure (e.g., h ¼ ðq, �Þ for a Newtonian fluid).Recall, that the shallow water equations only allow closure relations of the form rðh, u; hÞ:

Stiffness of the source term
The presence of friction implies stiffness of the system of equations that requires special attention in the construction of the numerical solver.The stiffness is associated with the relaxation time toward an equilibrium solution.Following Shi and Levermore (1996), the relaxation time T f is determined by the inverse of the largest absolute eigenvalue of the source term where k are the eigenvalues of the Jacobian @R @Q : For the shallow water system this is a known and well-investigated issue (see e.g., Shi and Levermore (1996); Xilin and Liang (2018)).For a Chezy type friction model one obtains a friction relaxation time scale Carrying over a similar analysis to the shallow moment system allows us to analyze whether the stiffness of the shallow moment hierarchy changes as we turn toward higher levels.
Assuming a Chezy-type bottom and Newtonian bulk friction closure with e x ¼ 0 yields the following eigenvalues of the Jacobian of the source term for systems up to level 6: We observe a cascaded structure, in which the lower-level eigenvalues are included in Equation ( 24) by taking the first l þ 2 rows of the vector, where l represents the chosen level.
We furthermore find that the eigenvalues grow by � 1 h in the bottom friction term and by � 1 h 2 for the bulk friction.For vanishing water depths, the Newtonian bulk friction term will therefore outgrow the contribution of the bottom friction quickly and the system will become stiff regardless of the values a i .
In addition to that, the eigenvalues scale with larger valued coefficients for an increasing shallow moment level.While the coefficients of a i might cancel each other, a pessimistic estimation (e.g., a i > 0 8i) also contributes to faster relaxation times of the friction force and therefore yields an even stiffer system.
Overall, we can conclude that the shallow moment system can potentially become very stiff for higher levels and small water depths, even compared to the already considered stiff shallow water equations.Therefore numerical algorithms for solving stiff equations must be taken into account when solving shallow moment systems.

Comparison to original VAM model
Here, we briefly want to discuss the similarities and differences of the shallow moment method with the VAM approach as presented by Ghamry (2000).Both approaches have in common that they are derived from the incompressible Navier-Stokes equations and use expansions of the velocities in order to compute moments of the mass and momentum balance in order to derive evolution equations.In the limit, using a constant expansion of the horizontal velocities (u, v), both models reduce to the shallow water equations (with different source terms).
While the shallow moment model makes assumptions on the flow conditions a priori by performing a dimensional analysis and neglecting terms of higher order, the VAM model starts with the non-reduced equations and makes assumptions on the flow conditions based on commonly used approximations, like the Boussinesq model which neglects laminar stresses.The shallow moment model on the other hand did not include normal deviatoric shear stresses based on arguments during the dimensional analysis, resulting in a purely hydrostatic pressure model and with a basal shear stresses r xz and r yz .
The expansion of the horizontal velocities in the VAM model takes the form where f 1 , f 2 , g 2 , g 2 are constraint by the relations The first constraint enforced that (u 0 , v 0 ) is indeed a mean velocity and the second constraint requires that (u 1 , v 1 ) introduce no net transverse transport.It is now interesting to observe that the Legendre polynomials considered for the shallow moment method fall into this class of ansatz functions if we define f 1 ¼ f 2 ¼ / 0 and g 1 ¼ g 2 ¼ / k for k � 1: This is true since / 0 ¼ 1 and the set f/ k , k ¼ 0, :::, Ng describes an orthogonal basis.
Theoretically, both, the VAM and the shallow moment approach therefore allow for a hierarchical development of evolution equations, since they both consider arbitrarily accurate representations of the vertical velocities u, v using a suitable basis.However, the VAM approach relies on a derivation of the moments by 'by hand', since the PDE is not transformed into reference coordinates and therefore the presented models are only available for linear and quadratic velocities u, v.The shallow moment approach is built upon the idea to 'automatically' generate the hierarchical evolution equations by using moment projections.It would be interesting to study if a similar approach could be taken for the VAM model.
The final system structure that needs to be solved is very simple in case of the shallow moment method, as already displayed in Equation (17).
The system for the VAM model takes the form The difference is that the time derivative acts upon a non-linear function wðQÞ: This makes solving the system more complex.
Additionally, if one compares both models with a linear expansion in the velocities u, v, the shallow moment method would result in 5 evolution equations (if we neglect the trivial equation for the bottom topography since it is not strictly required), compared to 10 equations in the VAM model.
Overall, the VAM model can therefore be considered as a model that is less restrictive in terms of the considered contributions of the stress tensor, resulting in a more complex model.The key difference of the shallow moment approach is its assumption of a hydrostatic pressure relation and it therefore neglects other normal stresses.This results in a significantly simpler set of equations.Since such moment models are typically compared to the shallow water equations, it is important that the final model is computationally cheap while containing the physically relevant information not available in lower-order systems.The shallow moment system can therefore be considered as a computationally very efficient choice for situations where non-hydrostatic effects can be neglected while still providing an adaptable degree of information about the vertical velocity structure.

Numerics
Is the goal of this section to present a simple and robust first-order finite volume scheme to solve system (17).
In the literature on finite volume methods, there already exists a vast selection of methods designed to overcome certain challenges, like well-balanced solvers, able to represent steady states; positivity preserving methods and solvers able to handle non-conservative products or stiff source terms.However, a lot of these methods are derived for the one-dimensional setting only or address a subset of the above requirements.The presented method will be able to handle stiff source terms and well-balanced with respect to the lake-at-rest solution.The method is written in a way applicable to unstructured two-dimensional meshes and can easily be implemented from an existing code base for hyperbolic PDEs without nonconservative contributions.
Our numerical approach is based on a first-order fractional step method by splitting the complete PDE into a PDE without source term and an ordinary differential equation (ODE) for solving said source term (Leveque and Yee 1990;Toro 2009;Holden et al. 2010).While the PDE is then solved explicitly using a first-order one-step method, the ODE can be solved in an implicit manner (or using higher-order explicit ODE integrators) to handle stiff source terms.The finite volume scheme implemented is able to handle non-conservative terms based on a class of path-conservative schemes using a weak formulation of the PDE introduced by Toumi (1992).The relevant literature for our particular implementation is based upon (Castro et al. 2006(Castro et al. , 2009a(Castro et al. , 2009b;;Canestrelli et al. 2009).

Treatment of stiff source terms
The partial differential equation ( 17) with source term and initial condition Q 0 reads Splitting Equation ( 26), e.g., following (Toro 2009), we obtain a PDE for the balance laws without the stiff source term and an ODE for the source term where Q denotes the solution of (27) after each iteration n.
The ODE ( 27) can now be solved, e.g., by using an implicit Euler scheme.Notice that solving the ODE implicitly for time step n with time step size dt only requires local information Qi and can therefore be performed independently for each element i.The linear system reads with the Jacobian of the source term denoted as @R @ Q ð Qi Þ: The size of the linear system is ð2N þ 2Þ � ð2N þ 2Þ: Since we are only interested in a first-order scheme, iteratively solving first for (27) followed by solving (28) is sufficient to compute Q nþ1 : An extension to higher order is possible, e.g., by using a Strang splitting approach (Leveque 2012).

Finite volume scheme
In order to derive our first-order finite volume scheme to solve (27), we now quickly recall the general derivation for hyperbolic systems with non-conservative products for unstructured meshes, using the path-conservative formulation proposed by (Castro et al. 2006(Castro et al. , 2009a(Castro et al. , 2009b)).
By rewriting the PDE ( 27) in its quasi-linear form by expanding @ x F y ðQÞ ¼ @F y @Q @ x Q and @ y F y ðQÞ ¼ where the quasilinear matrix is defined as A multi-dimensional one-step finite volume scheme can then be written as where jV i j denotes the area of the element i and jE ij j denotes the edge length of edge j of element i. n edges denotes the total number of edges per element.
where we used the segment path In order to obtain a Rusanov-type non-conservative scheme, we choose the dissipation term r ij to be r ij ¼ jkj max I WB , where jkj max is the maximal absolute eigenvalue and I WB is the modified identify matrix 1 0 ::: 0 1 0 1 ::: 0 0 . . . . . . . . .0 0 ::: 1 0 0 0 ::: 0 0 Setting the last row to zero ensures that the trivial equation for the bedload ( 5) is not violated.Combined with the additional one in the first row the method is also well-balanced with respect to the lake-at-rest solution.
For the practical implementation, we want to rewrite the current update scheme (32) to recover a form containing the classical flux function separated from contributions of the non-conservative terms.This can be done by manipulating expression (33) by using the definition of the quasi-linear matrix (31).For convenience, we also like to introduce the definition F :¼ ðF x , F y Þ T : Now, we can obtain the formulation based on Stokes' theorem, we can finally recast our update scheme (32) into the form In the above formulation (35), the first three contributions are identical to the implementation of a classical Rusanov scheme without the presence of non-conservative terms.Only the last term contains the non-conservative contribution.The term is computed using a Gaussian quadrature of third order in our implementation.

Treatment of boundary conditions
For the treatment of boundary conditions, we are following standard ideas (Leveque 2012;Toro 2009) to implement reflective, transmissive, and inflow conditions using a ghost-cell approach and in analogy to the boundary conditions used for the shallow water equations.A rigorous boundary condition analysis in order to define a well-posed problem is not yet available.However, the boundary treatment is in agreement with the boundary conditions proposed in the Vertically Averaged Moment (VAM) model proposed by Ghamry and Steffler (2002), an alternative model containing linear expansions of the horizontal velocities as well as dispersive terms.
The wall boundary or reflective boundary is implemented as a fixed, reflective and impermeable boundary.The data from the interior of the domain is mirrored to the ghost cell and the normal velocity components change sign.In case of the shallow moment model, velocity components refer to all projections of the velocities ðũ, ṽÞ The ghost cells are defined by for k ¼ 0, :::, N: The transmissive boundary condition can be implemented by mirroring the interior cell values to the ghost cell The inflow and outflow boundary conditions describe certain fields by providing data for the ghost cells.For a subcritical flow regime, we want to specify the mass influx where d a k and d b k are functions prescribing the boundary values.All remaining fields are implemented using a transmissive boundary.

Steffler's experiment
In this section, we apply the shallow moment model to simulate Steffler's experiment.Steffler's experiment conducted in 1984(Steffler 1984) observes the flow of water through an open, curved channel.The flow through a curved channel is particularly interesting for the study of the shallow moment system because it includes effects like superelevation and secondary flow.The latter is not representable by the shallow water system.Additionally, Steffler measured the vertical velocity profile of the flow, information that is not commonly available in geophysical experimental setups.Lastly, Steffler also simulated the experiment with a depth-averaged, dispersive model using linear reconstructions of the horizontal velocities (Ghamry andSteffler 2002, 2005;Vasquez, Millar, and Steffler 2011), already introduced as the VAM model.Both, the experimental data as well as the VAM model provide a valuable basis for comparisons with the shallow moment system.
The experimental setup is depicted in Figure 2, showing a 270 � curved rectangular flume with moderate curvature of R c =2b ¼ 3:4: The flume had a width of 2b ¼ 1:07m and the radius of curvature to the centerline of the channel was R c ¼ 3:66 m.The bed slope was kept constant at a value of 0.00083.The inflow of water starts at the bottom left into a 6 m long straight inflow section.After the 270 � bend, the water flows into a dropout after a 2.5 m long straight outflow section.The flow condition for the run used in our comparison was conducted with a velocity of 0.36 m/s and a depth of 0.061 m, resulting in a subcritical flow with Froude number 0.491.The boundary conditions are chosen similarly to Ghamry and Steffler (2002), namely prescribing the discharge at the inflow while setting all higher-order coefficients to zero.At the outflow, the height was enforced to be h ¼ 0:061m: The wall boundaries are implemented as described in Section 5.3.All necessary settings, in particular the values for all parameters, are presented in Table 1.
In Figure 2, we marked two different cross sections, located at a 90 � and 270 � angle of the curved channel.Those sections are used as reference  points for our following discussion.In each of these cross sections, we chose two sample points, one toward each side of the channel walls, located at R60:4b:

Secondary flow
By reconstructing the vertical velocity w using relation ( 6), we can compute streamlines for the two cross sections mentioned above.As already hinted at before, the circulation perpendicular to the main flow direction (longitudinal direction) is called secondary flow and cannot be represented using a shallow water model, due to the uniformity of the vertical velocity profile, visible in Figure 3.
For the level 1 shallow moment system, the secondary flow effect is clearly visible in the cross-sections of the flume.The rotation is always such that the water flows to the outside wall of the channel at the surface of the flow.A physical explanation for this behavior could be that the water is pushed outwards due to the centrifugal forces induced by the curved geometry.Some part of this generated momentum is converted into a rise of the water surface height at the outer wall, an effect called super-elevation, counteracted by gravity.The bottom friction lessens the momentum generated in the lower heights of the flow, such that parts of the momentum of the top layers can overcome the momentum in the lower layers resulting in a circulation direction dictated by the velocity at the surface.
Note that in some streamline figures shown above, there appears to be a flow through the wall boundaries.This, however, is not the case, as the normal velocity close to the wall approaches zero, as can be seen from the magnitude of the cross-section velocities depicted in color.
Since Figure 3 only gives a qualitative impression of the transversal flow profiles, the profiles of the transversal velocity in x-direction ũ and z-direction w are depicted in Figure 4.Note that the velocity profiles are plotted for a position in the middle of the cross-section, where the velocities w in z-direction are two orders of magnitude smaller than the horizontal velocities ũ: For both cross-sections located at 90 � and 270 � , the velocity profiles give expected results in the sense that for an increasing number of levels, the profiles seem to converge toward a specific profile.While the shapes of the first three levels are still clearly distinguishable, the profiles look very similar starting from level 4 onward.The velocity profile in z-direction shows the unexpected behavior that the deviation of the level 1 system compared � cross-section respectively.The rows show the velocity profiles u(z) and w(z) in x and z direction.
to the level 6 system is greatest.This indicates that the error in the moment hierarchy does not necessarily decrease monotonically, at least for systems with a low number of moments.A possible explanation for this behavior can be seen by looking at the definition of the velocity wðfÞ given in ( 6).This expression contains the term which is zero for a constant velocity profile at level 0. For level 1, where the gradient of the velocity profile is constant, the contribution is non-zero and potentially larger than the result for higher moments, where the vertical velocity profiles ðũ, ṽÞ become almost constant for f > 0:4 (see Figures 4 and 7).

Superelevation
In Figure 5, the surface heights are plotted over the cross-section at the two positions.The increase of the surface height toward the outer wall demonstrates the effect of superelevation.Comparing the 90 � and 270 � cross-section shows that the gradient of the water surface steepens throughout the flow in the curved section of the channel.While the 270 � cross-section qualitatively converged with level 2, it needs a system with level four to gain a similar accuracy at the 90 � cross-section.

Redistribution of the velocity profile
To investigate different shapes and therefore the redistribution of the vertical velocity profile at different locations, we visualized the profiles for the two cross-sections at 90 � and 270 � close to the inside and outside wall respectively.
In Figure 6, the transversal velocity profile u(z) is shown while Figure 7 displays the longitudinal direction v(z).At first glance, the transversal velocity profile shows very similar profiles for level 4 and higher, while the lower levels already give a reasonable approximation of the velocity profile and the mean value seems to be accurately represented.
The level 1 system tends to overestimate the velocity in the upper parts of the flow height.Interestingly, the inaccuracy of the level 1 and level 2 system is more prominent close to the outer wall.Here, it is especially obvious that in some situations a linear velocity profile is not able to capture the relevant vertical information and rather drastically overestimates the velocities at the upper section of the flow.
Having a closer look at the actual shape of the profile for higher moment systems at the different locations, it can be observed that the shape differs � and 270 � cross-section respectively.The rows show the transversal velocity profiles u inner wall ðzÞ and u outer wall ðzÞ: mostly in the lower section of the flow profile f < 0:3, while all profiles of higher level approach a roughly constant shape close to the surface f ! 1, due to the no-stress boundary condition.Since it is the shape at the bottom of the flow that is interesting, e.g., for computing a gradient to approximate shear stresses, the seemingly small changes in the lower section may become significant.
The longitudinal velocities show a more interesting behavior, since low moment approximations fail to accurately describe the mean velocity (compared to the mean of the highest moment) in some cases.Especially noteworthy is the fact that for the longitudinal velocity at the inner wall at 270 � , the mean velocity of the shallow water systems is more accurate than the mean velocity of the level 1 or level 2 system.We have no good explanation of why the moment hierarchy does not decrease in error for increasing levels in such low-moment systems.However, it is interesting to observe that this behavior also occurs at the inner wall for the cross-section at 90 � in a less pronounced fashion but it does not occur at the locations at the outer wall for both cross-sections.This again indicates the need for higher moment models to accurately capture important flow characteristics like the mean velocity.

Runtime comparison
In this section, we want to analyze the computational cost for an increasing number of moments (levels) in the shallow moment hierarchy.The data for the comparison is displayed in Table 2. Since the current implementation is an ad-hoc implementation of the numerics shown above without any optimization, the absolute results are normalized with respect to the shallow water simulation.The runtime increases by roughly 50% from one level to the next.One exception is level 4, where the runtime increased by ðþ84%Þ: We attribute this to the significantly higher than normal number of needed iterations until convergence.If this is an indication for a more expressive model needs to be investigated in the future.
The increase in runtime can be split up into two contributions, the first one being the increase of the runtime per iteration due to the larger system that needs to be solved.The second contribution is the number of iterations needed until convergence.For test cases where unsteady simulation results are important, a reasonable criteria for comparing the runtime of different levels of the moment hierarchy would be the norm: avg: runtime per unit time ¼ norm: avg: runtime per iteration norm: avg: time step size , which increases slower han the normalized total runtime due to the influence of the number of iterations needed to reach a steady solution.The increase of the normalized average runtime per unit time is mostly below 30% from one level to the next ðþ46% : Surprisingly, the average time step size, which scales with the largest absolute eigenvalue of the system in each iteration, starts to increase again after a first drop from level 0 to level 1.
In the current implementation, all eigenvalues are computed numerically, taking up a large portion of the overall computational cost.The construction of estimates for the eigenvalues based, e.g., on the level 1 system, where eigenvalues can be computed analytically, may provide a significant speed-up.Another big part of the computational cost stems from implicitly solving the source term.Changing the numerics to more modern approaches for handling stiff source terms, e.g., following (Dumbser, Enaux, and Toro 2008) will also lower the total computational cost.
While the shallow moment method has an increased computational cost for an increasing level of the hierarchy, the overall cost does not explode unreasonably.Therefore low-order systems will probably still be in an acceptable cost range for many applications.Since the cost of increasing the hierarchy by one level stays roughly constant, it would be very interesting to analyze if the error of the model compared to a fully three-dimensional simulation would also decrease at a similar rate and if there is an optimal level (in terms of computational cost versus accuracy) for certain applications.This, however, will be a topic of future research.

Comparison with reference data
The previous discussion focused on how the shallow moment hierarchy can be used to determine higher order effects that are typically not resolvable in the shallow water limit (N ¼ 0).
In this section, we will conduct a preliminary direct comparison of our numerical results with selected experimental and numerical data available from Steffler's experiment in Ghamry and Steffler (2005).In Figure 8, the height profiles of the cross-section at the 90 � and 270 � bend of the channel are compared with the experimental results.Note that a direct quantitative comparison suffers from major uncertainties in the resolution of the digitized experimental data from the paper (we could not access the data by any other means).The qualitative trend nonetheless shows that higher level shallow moment equations result in a better fit of the experimental results.The level 1 shallow moment model performs well for the 270 � cross-section while it does not very accurately match the data for the 90 � cross-section.The next Figure 9 shows a comparison of the longitudinal and transversal velocity profile at the 270 � bend at the outer wall of the channel.The plot shows a comparison with the experimental data as well as with the numerical solution obtained from the VAM model.Unlike in the previous comparison of the flow heights, the higher level shallow moment system does not provide a better fit to the given data.The level 1 system provides the best fit out of all shallow moment simulations.The higher-level simulations generate profiles that more accurately satisfy the free-surface boundary condition while the experimental data is not well captured.
The authors attribute this mismatch to the following factors: � The longitudinal velocity profile shows that the shallow water (level 0) and level 6 system do not capture the correct mean of the experimentally measured velocity.Shifting the velocity profile to the measured mean velocity would already greatly improve the result.The cause of this behavior is likely to be attributed to a non-optimal calibration of the model.This is something, however, that cannot easily be accounted for given the data source available to us, and has to be looked at more carefully in future work.(2005).The comparison is made at the 270 � bend at the outer section of the channel.
The experimental and numerical results from Ghamry and Steffler (2005) were digitalized from the publication figure which may negatively affect the resolution of the data.
� The Chezy coefficient was chosen based on the configuration of the simulation results of Ghamry and Steffler (2005).In our approach, however, we additionally provided a bulk friction term associated with a Newtonian friction model for water.Since this bulk friction term is not present in classical shallow water models, the calibration of the Chezy coefficient implicitly takes the resistance of the water body into account.Therefore our model would currently account for the resistance of the water body in both terms.This again implies a more careful look into how these types of models are empirically closed via friction relations and how they will be calibrated against measured data.One idea for the future would be to treat the Newtonian friction coefficient as an a priori measurable physical parameter in the system.As a result, the bottom friction closure becomes the only heuristical component of the bottom friction relation that needs to be calibrated based on observations.� The current implementation relies on a weak coupling of the free-surface and bottom-friction boundary conditions.Alternatively, a strong coupling of the boundary conditions could be realized.This would allow for alternative models, e.g., using a no-slip boundary condition at the bottom instead of the Chezy friction model.This requires a careful investigation by itself and is outside the scope of the current work.It will be addressed in a future study.

Conclusion
In this article, we derived the shallow moment method, implemented it in two space dimensions, and showed a simulation result based on a curved channel flow following Steffler's experiment.We demonstrated that even the simplest model, namely the level-1-system, can reproduce the secondary flow effect in the cross-section of the curved flume, an effect not visible in a shallow-water framework.We also compared the influence of the shallow moments hierarchies level on the resulting water surface elevation and the superelevation effect in the curved section of the flume.Superelevation is visible regardless of the level, however, the absolute value of the surface height differs significantly for the level 0 and level 1 system from the higher-order results.This is directly linked to the large deviation of the longitudinal velocity profile for those two systems compared to the higher-order solutions since total discharge is fixed by the inflow boundary condition.
The redistribution of the velocity profile's shape at different positions in the channel, be it in the same cross-section or different locations entirely, is an observation already demonstrated in previous papers and the main reason why the shallow moment model was developed in the first place.This contribution again confirms the need for higher levels to accurately describe the changing character of the velocity profiles, even in the steady state where temporal changes are not considered.
The vertical velocity profiles tend to converge quickly toward a common shape by increasing the level.This already happens for a low number of levels.This was shown by a qualitative investigation of the velocity profiles, where the profiles for systems with a level higher than four only changed marginally, compared to the more drastic deviations in lower levels.
Another, maybe more surprising observation is that low-order systems, like the shallow water or level 1 system, are not able to accurately capture the mean of the longitudinal velocity profile.This observation is of course not new and the introduction of a shape-factor in the shallow water equations is a classical and simple approach to correct for this deviation (Hutter, Wang, and Pudasaini 2005).However, we would like to emphasize that the shape factor is typically constant and therefore not able to represent the transient or spatial changes of the velocity profiles.If the shallowmoment method can be used as an intermediary to construct a transient and space-dependent shape factor for the shallow-water equations remains one interesting idea for future research.
A comparison of the computational cost of the shallow moment method shows that each additional level increases the cost by roughly 30% for unsteady simulations.For steady-state simulations, the total cost will additionally increase since it will take longer to converge to a steady solution.
The current article focused on a qualitative assessment of the benefits of using a shallow moment approach compared to classical shallow water models and the variances between different levels of the moment hierarchy.A direct comparison with the experimental and numerical data documented in Ghamry and Steffler (2005) showed that it is challenging to improve the predictive accuracy in models that also contain empirical parameters that need to be calibrated based on observations.This is due to the fact that empirical closures used in classical process models tend to compensate for shortcomings that result from a lower fidelity level of the system.In order to fully benefit from an increased accuracy of the shallow moment system compared to the shallow water model we will therefore in the future need to disentangle a gain in predictive quality that stems from the model formulation itself and the contribution to model error that stems from the heuristics in the system along with its calibration quality.Proper calibration of the shallow moment method requires more effort than a correct estimation of the friction parameter as it can be done for the shallow water limit.Additional steps involve the selection of an appropriate projection of the friction model and its interplay with the availability of bulk friction models, e.g., for Newtonian fluids.While these mechanisms allow for richer model constructions, they require further research to establish robust models for different physical settings.
In future work, we would like to strengthen the quantitative assessment of the shallow moment system for real-world and laboratory problems like Steffler's experiment.The availability of a fully three-dimensional numerical reference solution would greatly ease the process of finding appropriate closure models for a rigorous assessment of the capability of the shallow moment cascade, which is currently not available.

Figure 1 .
Figure 1.Schematical depiction of the different flow mechanisms.

Figure 2 .
Figure 2. Setup of the numerical experiment designed by Steffler in 1984.The mesh used for the simulation is used as depicted in the figure with 2150 quadrilateral elements.

Figure 3 .
Figure 3. Streamlines for level 0 to 4 for the two cross-sections at 90 � and 270 � angle of the curved section of the flume.u,w indicate the velocities in x and z-direction.

Figure 4 .
Figure 4. Velocity profiles at the middle of the cross-sections.The columns show the 90 � and 270� cross-section respectively.The rows show the velocity profiles u(z) and w(z) in x and z direction.

Figure 5 .
Figure 5. Plots of the surface height elevation at two cross-sections for various levels of the moment hierarchy.Both cross sections show the super-elevation effect due to the curvature of the channel.

Figure 6 .
Figure 6.Redistribution of the vertical velocity profiles at the inner and outer wall (R60:4b) of the cross-sections.The columns show the 90� and 270 � cross-section respectively.The rows show the transversal velocity profiles u inner wall ðzÞ and u outer wall ðzÞ:

Figure 7 .
Figure 7. Longitudinal velocity profiles at the inner and outer wall (R60:4b) of the cross-sections.The columns show the 90� and 270 � cross-section respectively.The rows show the longitudinal velocity profiles v inner wall ðzÞ and v outer wall ðzÞ:

Figure 8 .
Figure 8.Comparison of the water height obtained from the shallow moment model with the experimental measurements from Steffler's experiment Ghamry and Steffler (2005).The experimental results were digitalized from the publication figure and the resolution of the measurements should only be considered for a qualitative comparison due to an insufficient resolution of the original graphic.

Figure 9 .
Figure 9.Comparison of the longitudinal velocity profile (left) and transversal velocity profile (right) obtained from the simulations of the shallow moment model and the experimental measurements and numerical results (VAM) from Steffler's experimentGhamry and Steffler (2005).The comparison is made at the 270 � bend at the outer section of the channel.

Table 1 .
Numerical test case configuration.

Table 2 .
Comparison of different runtime characteristics for different levels of the moment hierarchy.