A High‑Order Conservative Semi‑Lagrangian Solver for 3D Free Surface Flows with Sediment Transport on Voronoi Meshes

In this paper, we present a conservative semi-Lagrangian scheme designed for the numerical solution of 3D hydrostatic free surface flows involving sediment transport on unstructured Voronoi meshes. A high-order reconstruction procedure is employed for obtaining a piecewise polynomial representation of the velocity field and sediment concentration within each control volume. This is subsequently exploited for the numerical integration of the Lagrangian trajectories needed for the discretization of the nonlinear convective and viscous terms. The presented method is fully conservative by construction, since the transported quantity or the vector field is integrated for each cell over the deformed volume obtained at the foot of the characteristics that arises from all the vertexes defining the computational element. The semi-Lagrangian approach allows the numerical scheme to be unconditionally stable for what concerns the advection part of the governing equations. Furthermore, a semi-implicit discretization permits to relax the time step restriction due to the acoustic impedance, hence yielding a stability condition which depends only on the explicit discretization of the viscous terms. A decoupled approach is then employed for the hydrostatic fluid solver and the transport of suspended sediment, which is assumed to be passive. The accuracy and the robustness of the resulting conservative semi-Lagrangian scheme are assessed through a suite of test cases and compared against the analytical solution whenever is known. The new numerical scheme can reach up to fourth order of accuracy on general orthogonal meshes composed by Voronoi polygons.


Introduction
The transport of sediment in rivers, lakes, or shores is particularly relevant in hydraulic engineering for several different reasons [2]. Erosion and deposition zones are generated by sediments which typically cause the modification of the river beds or that of the shores over time. Sediments have an influence on energy production in dams, and they may cause damages and serious risks for the environment, and they are responsible for delivering nutrients and, as opposite, contaminants. For the above reasons, the modeling of such flows, as well as the understanding of their dynamics, is of paramount importance in applications [38]. From the mathematical modeling point of view, multi-phase approaches are typically considered [3,11,18,30]. They consist in coupling the Navier-Stokes (NS) incompressible hydrodynamic descriptions for the liquid phase with advection-diffusion equations for the concentration of sediments. Alternative modeling approaches treat the sediments as a fluid, which turns out to be a reasonable choice when high-density sediment concentrations coexist with the fluid phases [39]. If one is instead interested in the understanding of the microscopic modifications of the fluid, each individual sediment is likely to be modeled via a microscopic description [29].
In this work, we consider three-dimensional free surface flows (typically water at the interface with air) modeled by the incompressible NS equations with a hydrostatic pressure distribution coupled with a scalar equation for the sediment transport. Sediments are supposed to be passive media. The equation for the solid-phase particles moves with speed determined by the liquid phase and an anistropic diffusion is added to the system. We focus on sediments in suspension and, thus, the main hypothesis is that the concentration of the solid phase is quite low. Additionally, sediments travel much more slowly compared to the fluid, and hence, they are supposed not to modify or affect its motion. This, in particular, can be shown to be true in slow current conditions and it can be inferred through a characteristic analysis of the problem [2]. Such analysis shows also that the variation of the ground occurs with a speed which is at least one order of magnitude lower than the variation of the free-stream flow. For this reason, we do not consider a morphodynamic model in this work, i.e., we do not consider the situation in which sediments may modify the bottom over time. This choice is also legitimated by the fact that the solid phase is dilute compared to the fluid one.
The numerical scheme chosen for discretizing the liquid phase takes inspiration from a series of papers by Casulli and co-authors [19,20,23,24,43]. More in detail, we focus on a semi-implicit time approach for the solution of hydrostatic free surface flows. The main advantage of such the method is that it permits to get larger stability domains while maintaining a reasonable computational cost. In fact, one of the main challenges related to this problem is due to the time step restriction induced by the fast scale coming from the viscous terms in both the horizontal and the vertical directions, together with the free surface speed. They typically cause a very severe CFL stability condition which penalizes explicit numerical schemes and which should be avoided if possible. The pressure terms, the velocity, and the vertical viscosity are discretized implicitly, while the nonlinear convective and horizontal viscous terms are discretized explicitly. The computational domain is discretized on the horizontal x-y plane by an orthogonal unstructured staggered grid where the normal velocity components are defined at the cell interfaces, whereas the pressure is evaluated at the center of each control volumes. In the vertical direction, we employ a standard z-layer discretization. A finite-difference scheme is used for the discretization of the momentum equations, and a mass conservative finite volume method is used for the time evolution of the free surface elevation.
The method employs polygonal Voronoi meshes, and makes use of a polynomial-based highorder reconstruction, along the lines of [9,28,33]. Concerning the remaining explicit terms in the NS equations and the solid phase, a semi-Lagrangian method is adopted [25,32,36,40]. Semi-Lagrangian methods require the integration of the fluid trajectories backward in time. This can be done in several ways. Here, we adopt a conservative approach [27,46] and a suitable interpolation of the vector velocity field that is provided by the high-order reconstruction procedure. Semi-Lagrangian methods see their origin for the numerical weather prediction [5,6,37,41,44,45]. Nowadays, they can be found in environmental engineering applications, such as free surface flows in rivers and oceans [27,46] as well as in plasma physics [4,17,26], in applications to image processing [15,16] or for solving the Hamilton-Jacobi equations [14]. Conservative semi-Lagrangian methods are described for instance in [4,31,36,40]. A fully space-time strategy with a semi-Lagrangian scheme on unstructured meshes has been recently proposed in [8]. In addition to the inherent conservative property, slope limiters can be introduced in the reconstruction to ensure preservation of specific properties such as positivity and monotonicity. In this work, we focus on a high-order semi-Lagrangian scheme which ensures conservation, hence presenting several advantages. Specifically, both momentum and sediment concentration will be proven to be conserved upon numerical evidences.
The outline of this article is as follows: Sect. 2 presents the physical model with the governing equations considered in this work, while Sect. 3 contains the details of the mesh, the definition of the discrete quantities, and the time discretization for both the fluid and the solid phases. Section 4 is dedicated to the high-order polynomial reconstruction and on the details of the semi-Lagrangian approach. In Sect. 5, a wide set of benchmark test problems is run to assess the accuracy and the validity of the new scheme. Finally, conclusions and outlook to future work are given in Sect. 6.

Governing Equations
The fluid phase is described by the three-dimensional hydrostatic free surface equations [22]. They are derived starting from the incompressible NS equations in the additional hypothesis that the vertical acceleration as well as the vertical viscosity forces are small when compared to the gravity acceleration and to the pressure gradient in the vertical direction. In Cartesian coordinates, these equations read where Eq. (1) expresses the mass conservation and, due to the hypothesis t = 0 on which the model relies, the equivalent incompressibility constraint. Momentum conservation in the x-and y-directions is modeled by Eqs. (2) and (3), respectively. The last Eq. (4) gives the time evolution of the water surface elevation, (x, y, t) , measured from the zero reference level located at z = 0 , as shown in Fig. 1. This equation is obtained integrating the continuity Eq. (1) over the water depth with the kinematic conditions where, u s , v s , and w s represent the velocity components at the free surface, while u b , v b , and w b are the corresponding components at the bottom. This is identified by the function h(x, y), measured from the reference plane z = 0 (see Fig. 1 for a visual explanation). The condition imposed to the bottom states that the velocity component perpendicular to the solid boundaries must vanish, and hence, no physical flux is allowed to cross the bottom. The total water depth is denoted by H(x, y, t) = h(x, y) + (x, y, t) . The hydrostatic approximation leads, neglecting the viscous terms in the third momentum equation, to p z = −g , i.e., the derivative of the pressure field in the z-direction is matched by the gravity acceleration g = 9.81 . Consequently, the normalized pressure is given by p(x, y, z, t) = p a + g( − z) with p a the atmospheric pressure. Finally, is the kinematic viscosity coefficient which can be taken as constant, or conversely, it can be determined starting from a specific turbulence model.
Mass conservation of a scalar variable is described by means of the following relation [13]: where c(x, y, z, t) denotes the concentration of a scalar transported quantity: salinity, temperature, suspended sediment, or any passive substance which may be relevant. The quantity w s appearing in Eq. (7) is the settling velocity: it is assumed constant and it has to be different from zero in the sediment transport case. The turbulent diffusion coefficients along the horizontal and the vertical directions are denoted by e h and e v , respectively.

Numerical Scheme for the Fluid and the Solid Phases
We introduce now the scheme; taking inspiration from [24], we employ to resolve the fluid phase described by Eqs. (1)-(4) and we give the details of the discretization of Eq. (7) for the solid phase. For the fluid phase, we propose a semi-implicit method designed to work on an unstructured grid, while unknown quantities are defined with a staggered approach. Before giving the details of the schemes, in the next section, we briefly discuss the employed mesh.

Computational Grid: an Unstructured Staggered Voronoi Mesh Approach
Let us first focus on the Voronoi tessellation [9,10] for the meshing of the computational domain . The class of unstructured grid is very important in fluid dynamics, because it allows to reproduce complicated geometry as the channel bends. The final grid is obtained carrying out two different steps, namely (i) the generation of a grid lying on the x-y plane xy , composed by a Voronoi tessellation, and (ii) the extrusion of the horizontal mesh along the vertical direction to produce the three-dimensional domain. The choice of using a Voronoi mesh is motivated by the fact that it is orthogonal by construction: namely the segment connecting two centroids of two adjacent elements is orthogonal to the shared edge that they have in common. This provides an essential property for the finite-difference discretization of the momentum equations discussed next. Here, we limit us to construct conforming meshes composed by convex elements.
More precisely, the generation of the domain xy is composed by three different steps. We briefly recall them: • the starting point is a Delaunay triangulation, referred to as primary grid, that must contain triangles with angles no greater than 90 • of amplitude; • after specifying an arbitrary refinement factor , every triangle is split into a total number of ( +1)( +2) 2 sub-triangles, hence generating an isotropic refined Delaunay triangulation; • the final Voronoi tessellation is then obtained by linking the triangle circumcenters between them. Figure 2 shows the dual bond between primary triangulation and Voronoi grid, where the Voronoi centers correspond to the triangle vertexes, while the circumcenters coincide with the Voronoi nodes. Along the horizontal boundaries, we note an additional layer that arises from the construction of the dual Voronoi tessellation, whose thickness depends on the characteristic mesh spacing of the primary grid. As already mentioned, we furthermore require that all primary triangles are acute-angled to avoid an empty intersection between two centers. The final three-dimensional computational domain is made of right prisms whose bottom and top faces are Voronoi polygons and height is a vertical layer of thickness Δz , that is assumed to be the uniform distance between two consecutive horizontal domains xy along the z-direction. Figure 2 depicts the final computational grid obtained for a generic domain both on the x-y plane and for the vertical direction. Now, we refer to Fig. 3 to introduce the notation adopted over the x-y plane and along the vertical direction z. We call the single Voronoi polygon as cell or element, and we indicate it with P i ∈ N P , with N P being the number of total cells, and P i representing its boundary which is composed by a total number of N j edges denoted with j . For a generic side j , let L(j) and R(j) be its left and right elements, respectively, with their corresponding centers L(j) and R(j) that are connected by the straight line segment j . The left and right orientation is given by the sign function i,j , defined on every edge j : We discretize the vertical column water by a sequence of z-layers of thickness Δz n j,k or Δz n i,k , depending if we are referring to a generic cell P i or to an edge j , as depicted in Fig. 3. Generally, we adopt a uniform height which remains constant, apart from those layers which are crossed by either the bottom or free surface, that have to be adjusted in order to exactly match the vertical boundaries of the domain. Indeed, some layers can be turned off or on, according to the time evolution of the free surface. More precisely, a cell or an edge is said to be active only if they are wet. We indicate the layer containing the bottom with the index m j and that one crossed by the free surface with M n j , and hence, the total number of active layers on the generic edge j is given by N n z,j = M n j − m j + 1 . If no water is present on edge j , then N n z,j = 0 and the edge is dry. In this context, we consider solid transport due to suspended sediments, but erosion and deposition mechanisms are not taken into account, and thus, m j remains constant in time. Finally, in Fig. 4, we present the staggered definition of the variables appearing in the mathematical models (1)-(4) and (7): -the free surface n i is evaluated for every cell P i at center i ; -the bottom topography h j and the total water depth H n j are located on each edge j ; -instead of storing the horizontal velocity components u and v as two separate variables, they are replaced by the velocity u n j,k , defined at the intersection point between the segment j and the edge j over the x-y plane and at the midpoint of the kth layer along the vertical direction. It is, furthermore, defined in the normal direction with respect to the edge j , pointing in accordance with the grid orientation given by Eq. -the concentration c i,k lies on the center of polygon P i as the vertical velocity but at height z i,k .
We are now ready to give the details of the space and time discretization.

Space and Time Discretization
The space and time discretization for the fluid phase is given along the lines of [24] and reads as follows: In the above equations, velocity, pressure, and vertical viscosity are discretized implicitly, while the convective and diffusive terms are computed explicitly. One consequently expects that the stability of the scheme depends only on the discretization of the nonlinear convective and the horizontal viscosity terms. Moreover, as it can be seen from (9), the horizontal momentum Eqs. (2)-(3) are reduced to the sole solution of the velocity in the normal direction j for each computational edge j and each active layer k. A semi-implicit finite-difference discretization is adopted for the momentum, i.e., for u n+1 j,k , while a massconservative finite volume scheme is used for the evolution of the free surface n+1 i . In addition, | j | denotes the length of the jth edge, |P i | the area of the Voronoi polygon P i on the x-y plane, and Δt = t n+1 − t n the current time step. The term Fu n j,k indicates an explicit finite-difference operator for the discretization of the nonlinear convective and horizontal viscous components. The details of it are given in the next Sect. 4.2. Systems (9)-(10) are solved by inserting the momentum Eq. (9) into the free surface Eq. (10), yielding a linear system in which the pressure n+1 i is the only unknown. The matrix of this system is symmetric and positive definite, and thus, the conjugate gradient method can be used for efficiently computing the free surface elevation at the next time level. Once n+1 i has been determined for every i ∈ [1, N P ] , the horizontal normal velocities u n+1 j,k are updated from Eq. (9). The finite volume discretization of the continuity Eq. (1) is finally used to update the velocity in the vertical direction as hence ensuring mass conservation. This concludes the fluid-phase time evolution. We discuss now the discretization of the solid phase expressed by Eq. (7). This reads In the above equation, Fc n i,k is an explicit finite-difference operator for the discretization of the nonlinear convective and horizontal viscous terms analogous to the one used for determining the time evolution of the velocity field for the fluid phase. However, as opposite to (9), the unknown c n+1 i,k lies in the middle of each z-layer in the vertical direction and, on . the x-y plane, in the two-dimensional cell center, while the velocity field is defined on the edges j of the Voronoi elements on the same x-y plane.
Since the discrete operators Fu n j,k and Fc n i,k are given by a semi-Lagrangian approach [22], they do not require any particular CFL condition to be stable. Indeed, the time step limitations are only imposed by the diffusion coefficients for both phases in the horizontal direction. We then have where l min = min √ �P i � is the minimum incircle diameter of the entire computational mesh and (e h , ) are the horizontal diffusion coefficients of the solid and fluid phases, respectively. Let observe that Eqs. (9) and (10) are coupled. Thus, in principle, one needs to solve the resulting system by means of iterative methods. However, due to the semi-implicit structure imposed, the problem can be strongly simplified. Let us start by introducing the following notations: Using the above definitions, the equations for the time evolution of the fluid phase (9) and (10) can be rewritten in the compact matrix-vector form as where, for completeness, one has to introduce inside the matrix n j the bottom and free surface boundary conditions. These are, respectively, u j,m j − 1 Once the velocity and pressure fields at time n + 1 are known, the solution of the solid phase can be computed. Using the matrix notation as done for the fluid system, Eq. (12) becomes This corresponds once again to the solution of a tridiagonal linear system for each water column and Thomas algorithm is used to obtain the concentration distribution values c n+1 i,k at the next time level.

Interaction Between Fluid and Solid Phases
Equations (1)-(4) with (7) are only able to describe the evolution of free surface flows containing a scalar passive quantity advected at the speed of the fluid. This means that no interaction between fluid and solid phases has been considered so far apart from the introduction of the settling velocity w s which constitutes a first step towards the coupling of both fluid and sediments. The next step that we briefly detail in this paragraph is to introduce a term which causes the fluid to be affected by the sediments. In particular, we follow the lines of [12,34,35] where sediment-driven flows have been considered in a sort of shallow water model for the fluid phase. The main modification with respect to Eqs.
(1)-(4) consists in a new pressure gradient in the momentum equations which takes into account the presence of suspended sediments. The new model reads (14) n j n+1 j where R > 0 represents a constant coefficient that measures the interaction between water and sediment, while c is the concentration of sediments. Analytical expressions for R can be found in the literature, as, for instance, in [12,34,35], and typically, they are concerned with a function R = R( , s ) , which then depends on the water and sediment density, and s , respectively. Assuming the density of the water constant and fixed, R is monotone increasing with the density s . An explicit discretization of the correction term 1 + Rc in the horizontal momentum Eq. (9) yields The concentration is computed on the common edge j at the layer height z k by a simple arithmetic average between the neighbor elements R(j) and L(j). From this relation, it is easy to deduce that the greater the density of the sediments and the concentration are, the greater the effects of the correction term on the pressure are present. As a consequence, horizontal and vertical velocities are affected as well, thus involving a modification of the entire momentum field. From the numerical point of view, the pressure correction 1 + Rc can be conveniently embedded in the governing equations as a coefficient that multiplies the implicit pressure gradient given by n+1 R(j) − n+1 L(j) . In this way, the time pressure wave Eq. (16) can be solved following the same path detailed in the previous paragraph without further modifications. Finally, the interaction term can be neglected by simply assuming R = 0 , thus recovering the initial model with a passive sediment. Evidences about the effects of this correction term for the pressure in the modified momentum Eq. (20) will be numerically shown in Sect. 5.6.

High-Order Reconstruction Procedure
In this part, we describe how in practice the operators Fu n j,k and Fc n i,k , defined, respectively, in Eqs. (9) and (12), are computed. We recall that these terms discretize the convective terms (respectively, nonlinear and linear) and the horizontal viscosity terms present in the momentum Eqs. (2)-(3) and in the equation for the sediment transport (7). The idea is to use a semi-Lagrangian strategy for these parts of the system. This requires a backward integration of the Lagrangian trajectories to find the foot of the characteristic for each unknown and, at the same time, a reconstruction of the transported quantity to evaluate the unknown function at this position. The numerical scheme presented so far makes use of pointwise horizontal velocity components u n+1 j,k , vertical velocity components w n+1 i,k+ 1 2 , and concentrations c n+1 i,k . Therefore, in the next section, we discuss a reconstruction procedure which permits to obtain a high-order description of the velocity field and of the sediment concentration within each computational cell. The details of the semi-Lagrangian technique will be instead given in Sect. 4.2.

Three-Dimensional Polynomial Reconstruction and L2 Projection Operator
The high-order description of the velocity field and of the concentration distribution is obtained by means of a polynomial reconstruction procedure. Within each control active is thus defined a three-dimensional polynomial function of degree N, which allows to recover a reconstructed solution of accuracy (N + 1) at a generic point p = (x p , y p , z p ) , inside the cell where p is located. These polynomials are defined at each time step exploiting the knowledge of the solution at time n.
We start by discussing a three-dimensional decoupled reconstruction strategy introduced in [9,10]. This is based on two subsequent steps. The first step consists in reconstructing both the velocity field and the sediment concentration on the horizontal plane x-y for each active layer k producing the high-order two-dimensional polynomials . Successively, the second step is concerned with a reconstruction along the z-direction based on a Lagrange interpolation. This second procedure uses the two-dimensional polynomials obtained in the first step. In general, the choice of a decoupled reconstruction for the horizontal and the vertical flow directions allows the scheme to be computationally less expensive compared to a fully three-dimensional reconstruction, and for this reason, it is employed here.
The two-dimensional reconstruction polynomials in the x-y plane are expressed by a normalized Taylor series of degree N using the pointwise values furnished by the scheme at time n, that is, The term l i = √ �P i � denotes the normalizing factor necessary to avoid ill-conditioned reconstruction matrices (see [1] for details) and c = (x i , y i ) represents the centroid, that is the center of the Taylor expansion and coincides with the centroid i . Notice that the summation over all terms of the Taylor series explicitly writes and counts a total number of addends l = 1, ⋯ , N with N = (N+1)(N+2)

2
. These are nothing but the degrees of freedom of a polynomial of degree N along each spatial direction, and thus, it follows that the reconstruction polynomials (21)-(23) can be written in compact notation by adopting the following expansions: where Einstein summation convention is adopted implying summation over repeated indexes. The unknown expansion coefficients are given by the space derivatives of the polynomials to be constructed, In the above formulae, the basis functions are denoted by l and are piecewise polynomials belonging to ℙ N . Let observe that the degrees of freedom ̂ l i,k for the horizontal velocity count twice the number N , because two spatial directions are involved, according to (21). Let us note that, for the horizontal velocity and for the concentration, the polynomials are defined at the height level z n i,k , whereas the vertical velocity is defined at location z n i,k+ 1 2 (see Fig. 3 for the definition of the variables in the system).
The unknown coefficients (28)-(30) are then recovered by solving a linear system. This linear system is composed by the values of the functions to be reconstructed on a given stencil. For the horizontal normal velocity, the stencil S i is composed by a total number n S of edges j , while for the vertical velocity and for the concentration, we consider the entire cells P j to construct the stencil S P i . Typically, due to the unstructured nature of the mesh, the number of elements composing the stencil is bigger than the necessary minimum number N of unknown expansion coefficients required to obtain the formal order of accuracy (N + 1) . In particular, to perform our reconstruction, a safety factor of 2 is adopted, and hence, we set n S = 4N for Eq. (28) and n S = 2N for Eqs. (29)- (30), so that an overdetermined system is obtained. For those cells lying on physical boundaries, the reconstruction stencil results in a one-sided stencil composed by elements inside the computational domain. The conditions to be fulfilled by the reconstruction polynomials are that conservation holds true in each of the cells composing the stencils [7,28], where integrals are computed using Gaussian quadrature rules of suitable order of accuracy, see [42]. As already stated, by construction, the resulting system is overdetermined. Thus, a solution can be found, in general, by solving it in the least square sense. Unfortunately, this inevitably leads to loss of conservation issues, since Eqs. (31)-(33) cannot be anymore solved exactly. This problem can be overcome by adding a linear constraint to the system, The role of the above constraints is to require that the reconstructed polynomials exactly match at least the pointwise values of the unknowns on the polygon P i where they are defined. To enter the details, the unknown coefficients of the Taylor expansion to be determined can now be conveniently written in vector notation as This leads to three overdetermined linear systems which can be formulated in matrix form as where their relative linear constraints also appear in the second rows. In the above expressions, u , w , and c are the reconstruction matrices depending on the chosen Taylor expansion which can be pre-computed once for all offline (since they are simple geometric objects), whereas n u , n w , and n c are the right hand-side vectors containing the cell values of velocities and concentrations that change at every time step. The second rows in Eqs. (40)- (42) correspond to the matrix form of the linear constraint condition. Finally, using a constrained least square approach, one gets where , , and are the vectors of Lagrange multipliers used to enforce the linear constraints.
The reconstruction can be extended to the vertical direction as well, starting from the twodimensional polynomials on the horizontal plane and relying on one-dimensional Lagrange polynomials along the z-direction. Due to the staggered grid discretization, we have to deal with two different situations. In particular, we define two different reconstruction stencils consisting of (N + 1) values. The first, S i,k , is used for the horizontal velocities and the sediment concentration, and is defined at height z k , the second, S i,k+ 1 2 , is used for the vertical velocity and is defined at level z k+ 1 2 . Thus, we obtain a continuous three-dimensional representation of the solution as with where each stencil counts a total number of N + 1 layers corresponding to the number of degrees of freedom in one dimension ( d = 1).
The decoupled reconstruction procedure described above guarantees high computational efficiency, because it involves a more compact reconstruction stencil compared to a truly threedimensional reconstruction. Nevertheless, for evaluating the reconstructed value of a quantity (either velocity or concentration) at a generic point p , it would be much more convenient to deal with a polynomial directly defined for d = 3 , i.e., within each control volume V i,k . This is the reason why we introduce the following three-dimensional polynomials n, 3D i,k ( ) and The integrals appearing in the above formulae are again evaluated using Gaussian quadrature rules of suitable order of accuracy, while the values of the functions on the right-hand side of Eqs. (54)-(55) are computed at each Gauss point by means of the decoupled highorder reconstruction polynomials (46)-(47). The unknown expansion coefficients ̂ 3D,l i,k and ̂ 3D,l i,k can, therefore, be computed. We will show producing the convergence studies reported in Sect. 5 that the accuracy order of this new three-dimensional polynomials remains equal to (N + 1) , as in the case of polynomials obtained from the decoupled reconstruction.

The Conservative Semi-Lagrangian Operators
In this section, we detail the discrete operators Fu n j,k and Fc n i,k contained in Eqs. (9) and (12), which account for the convective and horizontal diffusive terms of the free surface NS model (2)-(3) and the sediment transport model (7), respectively. The semi-Lagrangian approach for these operators reads as To design a conservative scheme, the three-dimensional control volume V n i,k is transported backward in time by integrating the Lagrangian trajectories that depart from all the vertexes that define it, as shown in Fig. 5 for a simplified setting with d = 2.
In this way, the transported volume V L i,k becomes the control volume onto which the transported quantities have to be integrated, namely the velocity field and the sediment concentration. Specifically, in Eq. (56), the horizontal velocity is integrated along the edge L j of V L i,k at height z L i,k , while for the concentration (57), we consider the Voronoi polygon P L i defined again in V L i,k at level z L i,k . Furthermore, L j,k and L i,k denote, in the new control volumes V L i,k , the points where the velocity and the concentration are defined, respectively, at the intersection point between the segment j and the edge j over the x-y plane and at the center of the polygon P i in the starting volume. These quantities are readily known calculating the feet of the Lagrangian trajectories connected to the vertexes constituting the Voronoi volume and are defined by applying an appropriate numerical integration procedure to the transported control volumes. In particular, what is integrated are the continuous polynomial functions that define the velocity field and the concentration distribution at time n obtained from the three-dimensional reconstruction procedure presented in the previous Sect. 4.1, i.e., polynomials n,3D ( ) and n,3D ( ) . The other two terms appearing in Eqs. (56) and (57), ∇ 2 h and ∇ 2 h c , are suitable discretizations of the horizontal Laplacian operators regarding viscosity and diffusion at the end of the Lagrangian trajectory. For consistency, neglecting convection and horizontal diffusion, formulae (56) and (57) reduce to the identities Fu n j,k = u n j,k and Fc n i,k = c n i,k . Now, using the reconstructed velocity field, one can numerically integrate the Lagrangian trajectories connected to all the vertexes (m 1 , m 2 , ⋯) of the mesh element to compute  ( ; , Δ ) . In practice, the condition (Δ ;x) = is fixed in such a way that we look for the curve passing on the generic vertex of the control volume m 1 after the time step Δ . A graphical representation of the above situation is given in Fig. 5 for Δ = Δt . The system (58) is then approximated by means of a high-order Taylor expansion. This permits to write In the above expression, L corresponds then to the searched foot of the characteristic corresponding to the point . Here, to compute the Lagrangian trajectories, we make the hypothesis that ū is frozen at time . Extensions of the proposed method to the case of time reconstructions can be obtained following the path discussed in [8]. Under this hypothesis, then, at the aid of the chain rule, the time derivatives in Eq. (59) can be replaced by spatial derivatives that are readily computed from the high-order reconstruction polynomials and permit a direct computation of the foot of the characteristics. This gives in practice, and consequently, The spatial derivatives appearing in Eq. (60) can be easily evaluated relying on the highorder expansion (49). Indeed, the derivatives are entirely taken into account by the basis functions, that is, The order of accuracy for the Taylor expansion (59) is chosen to be the same of the order adopted in the reconstruction procedure, so that all high-order terms can be easily evaluated. The numerical integration stops either at the time = Δt or when the physical boundary of the computational domain has been crossed by the characteristic curve. In the latter case, as shown in Fig. 5, the shifted control volume is still computed and consequently subdivided into a physical and a ghost partition, which does not and does exceed the limits of the domain, respectively. The convective operators (56), (57) are formally evaluated in the same manner, integrating the prescribed boundary condition over the ghost partition and the reconstructed numerical solution over the physical partition. Specifically, the flow field is assigned in the case of Dirichlet boundary conditions, while, for periodic boundary conditions, the ghost cells are simply the cells located on the other side of the torus which, by definition, lies inside the domain. Finally, Neumann conditions require the computation of the derivatives, which are reconstructed using only elements inside the domain and no ghost cells are needed.
The time step Δt is given by the stability condition introduced with Eq. (13). In the case in which the sediment transport is considered, we assume ̄ ( ) = (u( ), v( ), w( ) − w s ) . It might happen that the shifted control volumes at the foot of the characteristics are degenerate or inverted, even though this is very rare for free surface applications. Nevertheless, a possible remedy to overcome this problem could be to introduce an a posteriori indicator which permits to detect which control volumes do not fulfill some regularity conditions. Then, one can recompute the solution in these volumes by reducing the order of the polynomial reconstruction. Alternatively, one can reduce the time step or remapping the solution over a new grid. The above ideas have not been explored in this work and they will be the subject of future investigations.
The diffusive components ∇ 2 h and ∇ 2 h c of Fu n j,k and Fc n i,k are defined at the end of the Lagrangian trajectory. For a generic scalar quantity , the discrete Laplacian is obtained first using Gauss' theorem, and successively relying on the exact solution of the generalized Riemann problem of the heat equation to discretize the integral Eq. (63). This yields where − k corresponds to the reconstructed value onto the boundary P i from within the element P i and + k is the reconstructed value again onto the boundary P i but from the neighbor of P i . This concludes the description of the numerical method for the fluid phase and the solid phase.

3 5 Numerical Results
In this section, we present several tests which are used to validate the numerical schemes illustrated previously. Specifically, we focus on analyzing the new conservative semi-Lagrangian schemes proposed to advect the NS equations and for computing the solution of the solid-phase equation. In details, we propose the following numerical tests.
1) Convergence test for the high-order polynomial reconstruction.
2) Convergence test for the semi-Lagrangian method in the zero diffusion case.
3) Convergence test for the semi-Lagrangian method in the zero advection case. 4) Stationary vortex flow for the NS equations. 5) Advection of a sphere of contaminant along a channel. 6) Introduction of sediment into a channel by a point source. 7) Sediment transport into a channel with settling velocity. 8) Gaussian pulse with water-sediment interaction.
If not stated otherwise, we always assume a flat bottom profile set at height z = 0 , so that the free surface elevation coincides with the total water depth.

Convergence Study
To check the correct implementation of the high-order 3D reconstruction for the sediment concentration presented in Sect. 4, we assume a smooth initial distribution that reads The computational domain is the unit cube = [0, 1] 3 that is discretized with a fixed number of active vertical layers N z = 50 . Instead along the x-y plane, each layer is paved with an unstructured Voronoi mesh and a sequence of consecutive values of a refinement factor are considered. We measure the L 1 , L 2 , and L ∞ error norms of the reconstruction as where s represents a generic variable for which one wants to measure the error, while s e ( ) and s h ( ) denote the exact and the numerical reconstructed solutions for a given mesh size h, respectively. The above integrals are computed by summing the contributions of all the Voronoi active cells V i,k ∈ where the index i ∈ [1, N P ] identifies the 2D Voronoi polygon and the index k ∈ [1, N z ] the vertical layer. Three-dimensional Gaussian quadrature formulae of appropriate order [42] are then applied. The numerical order of accuracy p is evaluated as follows: (65) c 0 (x, y, z) = sin 2 x cos 2 y cos 2 z . (66) The results are displayed in Table 1. They confirm that the expected theoretical accuracy order has been reached and they validate the high-order reconstruction procedure up to fourth order. A second convergence study has been implemented to test the effectiveness of the conservative semi-Lagrangian discretization proposed in this article. The initial concentration is the same as the one given in Eq. (65), while for the NS equations, we consider the following initial data: corresponding to a uniform flow in longitudinal direction with constant free surface elevation. In addition, we take the diffusion coefficients (e h , e v ) = 0 in the solid-phase equation. The exact solution of the NS equations for this set of initial data consists in a free stream with constant velocity , while the sediment concentration is advected at the speed of the fluid phase: c(x, y, z, t) = c 0 (x − ut, y, z) . For this problem, we measure the errors produced by the scheme for the solid phase and we report a convergence table analogous to the one presented previously for the polynomial reconstruction by fixing the time step to Δt = 0.1 .     Table 2 shows the obtained results for different polynomial degree N. The convergence rates match the theoretical prediction for all tested polynomial reconstructions. We now study the convergence for our numerical scheme for the solid phase in the case in which pure diffusion is present, i.e., (e h , e v ) ≠ 0 , while the fluid is at rest. In particular, we focus on the case in which the diffusion constant e h is different from zero in the longitudinal x-direction and zero otherwise. In this situation, our model for the sediment transport reduces to the one-dimensional heat equation for which we consider the exact solution: where g(x) is an arbitrary initial distribution for the concentration. In our case, we assign to g(x) the following Riemann initial datum:  where e h = 10 −4 . We consider a square domain in the x-y plane given by xy ∈ [−0.5, 0.5] 2 , while along the vertical direction, we have z ∈ [0, 1] . This domain is discretized with a uniform mesh spacing of Δz = 0.2 , equivalent at N z = 5 in the vertical direction, while different mesh sizes are considered on the plane x-y to produce the convergence results. Figure 6 depicts a one-dimensional cut of the solution at two different instant of time: namely at the beginning of the simulation for t = 10 and when the simulation ends at t = 200 . In both panels, the exact and the computed solutions are compared to each other. We clearly see a good matching between the two profiles even for long times. Figure 7 depicts the projection of the same numerical solution of Fig. 6 in different directions.
Finally, in Table 3, we report the norm of the errors and the convergence rates obtained by refining the mesh in the x-y plane. Again, the theoretical order of convergence is reached by the numerical scheme for all the tested reconstructions.

Stationary Vortex
In this part, we consider the NS model and we test the new semi-Lagrangian scheme in the nonlinear case which we use to advect the momentum equation. The test is a stationary vortex with a variable free surface produced by a balance between the pressure gradient and the centrifugal force [9], where viscous terms are neglected. Thanks to the radial symmetry of the problem, the momentum equation in polar coordinates is reduced to a balance between pressure forces and angular velocity. It can then be written as where r = √ x 2 + y 2 is the radius and u the angular velocity. Then, to construct an analytical solution, one needs to balance the centrifugal force. A possible choice consists in adopting the following initial distribution for the angular velocity: It is possible to find an explicit expression for the free surface replacing (73) into (72) and integrating over the radial direction. This leads to r . where we assume as integration constant the unperturbed free surface elevation sufficiently far from the vortex center, that is ∞ = 1 . Figure 8 plots this steady solution for the tangential velocity and for the free surface, respectively. Instead, Fig. 9b shows the same profile for the free surface in a 3D view. In this setting, the computational domain is a circle with r = 10 on the x-y plane, paved with N P = 9 511 polygonal elements, while along the vertical direction, we adopt N z = 10 vertical layers in the interval z ∈ [0, 1] . We can observe the 2D Voronoi mesh used for this test in Fig. 9a. The final time is t = 0.5 and the time step is chosen to be Δt = 0.1 . Figure 10 depicts the numerical results produced for the tangential velocity and the free surface at final time t = 0.5 . For all the three polynomial reconstructions that have been tested, a very good matching between exact and numerical solutions is achieved. Furthermore, a zoom of the free surface is proposed in Fig. 11 where we can observe that the polynomial reconstructions of high order give a better reproduction of the exact solution.

Advection of a Sphere of Contaminant Over a Channel
This   Figure 13 shows some cross-section views of the results at time t = 1.5 , while Fig. 14 at t = 3 . They are compared with the results produced by the numerical scheme presented in [10] in which the semi-Lagrangian method for the terms Fu n j,k and Fc n i,k was not conservative. The results underline how the conservative version of explicit terms better preserves the solution over the time evolution of the solution. This is of paramount importance when the long time behavior of the simulation is demanded.

Point Source Solution
Here, we consider the introduction of suspended sediments inside a flat rectangular channel by a point source. For the fluid phase, we suppose a uniform flow in the longitudinal direction, a hydrostatic pressure, and an inviscid fluid. The boundary conditions for the sediment equation correspond to the imposition of null flux across the channel boundaries where y L and y R represent the transverse coordinates of the left bank and the right bank, respectively, whereas b and 0 denote the constant elevation of bottom and initial free surface. The described problem has an analytical solution which reads with (76) e h c y

Sediment Transport with Settling Velocity
The last test problem is concerned with the transport of a sphere of sediment with a settling velocity w s ≠ 0 in (7). The computational domain is the channel = ∈ ℝ 3 |x ∈ [0, 4] ∧ y ∈ [0, 1] ∧ z ∈ [0, 1] , that is paved with a total number of N P = 3 176 Voronoi polygons and N z = 50 uniform layers of thickness Δz = 0.02 . The initial condition for the fluid phase is simply given by an inviscid fluid with the velocity 0 = (1, 0, 0) and the free surface elevation 0 = 1 , while the sediment concentration is assumed to be initially unitary and defined within a sphere of radius r s = 0.1 centered at s = (0.5, 0.5, 0.75) , by means of (75). Two different simulations are run up to the final time t = 3 with time step Δt = 0.1 , namely with e v = 0 and e v = 5 ⋅ 10 −3 , while e h = 0 in both cases. The settling velocity is w s = 0.25 , so that the sediment sphere reaches the bottom of the channel where concentration starts to grow, hence accumulating sediment. Figure 17 reports the result for both simulations at different output times, showing a slice along the x-z plane at y = 0.5 , thus following the transport and, eventually the diffusion, of the sediment. If diffusion is present, the concentration is not accumulating as much as in the pure advected configuration, and hence, the sediment remains suspended, as expected. Finally, Fig. 18 demonstrates that the sediment mass as well as the fluid mass is conserved by the novel semi-Lagrangian scheme.

Gaussian Pulse with Water-Sediment Interaction
In this last part, we validate the coupled water-sediment model introduced in Sect. 3.3. This model corresponds to case in which the modified momentum equations (20) defining the influence of sediments on the fluid flow are introduced. The setting of the problem is the following. The computational domain is = [−1, 1] 2 × [0, 1.5] , discretized with 30 z-layers along the vertical direction and paved with 7 377 Voronoi cells over the x-y plane where no-slip wall boundary conditions are imposed. The initial pressure distribution is depicted in Fig. 19a and is given by 0) , s = 1 , and A = 0.5 . The fluid is initially at rest and the horizontal viscous terms are neglected, while the vertical viscosity coefficient is chosen to be = 10 −1 . A sphere of sediments with radius r = 0.2 is located inside the fluid at (x c , y c , z c ) = (0.2, 0, 0.6) . A uniform initial concentration equal to c = 500 is assigned to the sphere and the settling velocity is w s = 0.25 . The initial concentration distribution is depicted in Fig. 19b through a cross-section of the entire fluid domain. We consider two different cases, namely R = 0 and R = 0.014 , in order to enhance the differences between passive and active sediment transportation. Both simulations are run up to the final time t fin = 0.5 . Results at different output times are shown to qualitatively highlight how the correction term acting on the pressure in the moment equations affects the fluid flow. Figures 20 and 21 compare the velocity flow field at different cross-sections on, respectively, the x-z and x-y planes. Figure 22 depicts different views on the horizontal plane for the free surface elevation, while corresponding one-dimensional plots of the free surface along the line y = 0 are shown in Fig. 23, giving evidences of the different behavior of the flow with R = 0 and R = 0.014 . In particular, one can observe that symmetry is lost due to the presence of sediments in the case of a coupled approach between the fluid phase and the solid phase, while it is maintained in the opposite situation by the numerical scheme. Finally, Fig. 24 reports the concentration distribution across the plane x-z at y = 0 . In this last figure, it is possible to notice how the sediments undergoes a slightly different distortion depending on the passive or active transport that has been modeled.

Conclusions
In this work, we have presented a new semi-Lagrangian method for solving the threedimensional free surface hydrostatic NS equations coupled with an advection-diffusion model for the sediment transport. The numerical scheme works on arbitrary polygonal meshes and its semi-implicit character permits to overcome severe stability conditions imposed by the problem under consideration. The new proposed method is fully conservative, since the sediments and the velocity field are integrated for each element where they are defined over the advected volume. The new deformed volume is obtained by computing the feet of the characteristics arising from the vertexes of each element. High-order reconstruction of both the sediment concentration and the velocity field allows the Lagrangian trajectories of the flow particles to be integrated, and the quantities at the feet of the characteristics to be properly interpolated. In a second part, we have shown that the theoretical order of accuracy is reached at numerical level by the scheme both for the fluid as well as for the solid phase and, at the same time, very good results in terms of precision are obtained. In particular, we have shown that the proposed semi-Lagrangian method outperforms other standard high-order semi-Lagrangian methods recently proposed in terms of precision and for long time simulations.
In the future, we would like to go deeper into the modeling aspects related to the transport of sediment realizing a fully morphodynamic model which also includes a proper treatment of wetting and drying fronts [21]. We want to take into account the influence of the sediments on the fluid phase due to deposition over the bottom and the non-negligible turbulence effects caused by the suspended particles over the liquid motion. Finally, further investigations are also planned to find a remedy for the appearance of degenerated transported volumes that can arise from the backward integration of characteristics in the semi-Lagrangian approach.