A New Mass-Conservative, Two-Dimensional, Semi-Implicit Numerical Scheme for the Solution of the Navier-Stokes Equations in Gravel Bed Rivers with Erodible Fine Sediments

: The selective trapping and erosion of ﬁne particles that occur in a gravel bed river have important consequences for its stream ecology, water quality, and overall sediment budgeting. This is particularly relevant in water bodies that experience periodic alternation between sediment supply-limited conditions and high sediment loads, such as downstream from a dam. While experimental efforts have been spent to investigate ﬁne sediment erosion and transport in gravel bed rivers, a comprehensive overview of the leading processes is hampered by the difﬁculties in performing ﬂow ﬁeld measurements below the gravel crest level. In this work, a new two-dimensional, semi-implicit numerical scheme for the solution of the Navier-Stokes equations in the presence of deposited and erodible sediment is presented, and tested against analytical solutions and performing numerical tests. The scheme is mass-conservative, computationally efﬁcient, and allows for a ﬁne discretization of the computational domain. Overall, this makes the model suitable to appreciate small-scales phenomena such as inter-grain circulation cells, thus offering a valid alternative to evaluate the shear stress distribution, on which erosion and transport processes depend, compared to traditional experimental approaches. In this work, we present proof-of-concept of the proposed model, while future research will focus on its extension to a three-dimensional and parallelized version, and on its application to real case studies.


Introduction
River damming produces alterations on the natural river functionality, both in the water discharge, as well as in the sediment transport and connectivity [1,2]. In particular, large dams [3] are estimated to trap more than 99 % of the sediments entering the reservoir [4]. This causes the progressive silting of the reservoir and inhibits the sediment load in the river flow downstream of the dam, thus altering the river morphology [5] and the aggradation/degradation dynamics that are closely linked to the balance between upstream sediment supply and local transport capacity conditions [6]. Sediment supply-limited conditions from upstream cause the (selective) erosion of finer particles from the granular bed, until the flow is unable to move the coarser grains and new equilibrium conditions are reached [7]. During this degradation process, the median size of the bed material progressively coarsens and the sediment transport rate decreases, leading to a process known as bed armoring [8][9][10]. Simulation (DNS) over rough bed configurations do exist (e.g., [33][34][35]), the inclusion of sediment active layers in fine-resolution hydrodynamics model is a relatively unexplored area of research.
In this study, we present and test a new semi-implicit numerical scheme for the solution of the two-dimensional Navier-Stokes equations, in which we included the possibility to easily simulate sediment entertainment and transport processes. The scheme, based on the method proposed by [36][37][38], is mass-conservative, computationally efficient, and able to solve the small-scale structures that characterize inter-grain flow field. In this study, we present proof-of-concept and preliminary results of this model as a first step towards its extension to a complete three-dimensional model coupled with a turbulence closure scheme. To this end, we focus on the validation of the proposed model against numerical tests and on showing its potential for applications in the context of fine sediment transport dynamics in gravel bed rivers.

Governing Equations
The governing equations are the two-dimensional incompressible Navier-Stokes, that are composed by the momentum conservation: ∂w ∂t and the incompressibility condition: ∂u ∂x where u, w are the components of the velocity field in the x− and z− direction, ν = µ/ρ is the kinematic viscosity coefficient, µ is the dynamic viscosity coefficient, ρ is the fluid density, p is the pressure term normalized by the fluid density, and g is the gravity acceleration. If the pressure is assumed to be locally hydrostatic, it can be expressed in terms of the atmospheric pressure p a and the hydraulic head η as: p = p a + g(η − z).
By means of Equation (3) it is possible to rewrite the momentum Equations (1) in terms of the hydraulic head η: ∂u ∂t ∂w ∂t Further introducing the mass conservation law for the suspended sediment concentration c (here defined as volume concentration, hence simplifying the mass conservation law assuming constant sediment density), and defining its relation with the sediment level S = S(x, z, t), two more equations need to be added to the system: ∂S where w c is the falling velocity of the suspended sediment; f Sc and f cS are two functions that regulate the mass transfer between the passive concentration c and the sediment level S. Note that in Equations (5) and (6) the sediment dynamics are intrinsically assumed to be driven only by the sedimentation/erosion process through the functions f Sc and f cS . The suspended concentration is expressed as volume of suspended sediment per unit of control volume (i.e., volume concentration relative to each computational cell), thus it can range in c ∈ [0, 1 − φ], where φ indicates the sediment porosity. Similarly, the sediment level S is defined per unit of volume, thus it can vary in the range S ∈ [0, 1], where S = 1 corresponds to a computational cell completely filled by the sediment. Note that the suspended sediment is assumed to be passively transported by the fluid, while the deposited sediment level actively interacts with the fluid dynamics as in general it may vary in time, thus changing the boundaries of the fluid domain. The exchange between suspended and deposited sediment must satisfy the mass conservation law:

Staggered Mesh
The Partial Differential Equations (PDE) system composed by Equations (2), (4) and (5) is numerically solved in a domain Ω(t) ⊆ Ω ∈ 2 discretized by a regular mesh Ω = i,k , for i = 1, . . . , I max and k = 1, . . . , K max . The main quantities are defined following a staggered approach. In particular, the hydraulic head η, the suspended concentration c and the sediment level S are defined in the cell centers, namely η = η ik , c = c ik , S = S ik . The horizontal velocity is defined in the center of the vertical edges, while the vertical velocity in the center of horizontal edges, namely, u = u i± 1 2 ,k and w = w i,k± 1 2 .
The system is solved for discrete times t = t 0 , t 1 , . . . , t n , . . ., hence the solution at a certain time t n is marked with the upper index n, namely ξ(t n ) := ξ n , while the time interval between two consecutive steps is ∆t n = t n+1 − t n . Note that the effective volume in each control cell can have different values due to the presence of the sediment, and can be computed as V i,k = ∆x∆z(1 − S ik ). Since S i,k ∈ [0, 1], the volume computed in this way satisfies 0 ≤ V i,k ≤ ∆x∆z. The effective edge length can change as well, being defined above the sediment top: for any time t n the effective vertical and horizontal edges lengths are defined as ∆z n i± 1 2 ,k and ∆x n i,k± 1 2 , see Figure 1. For non-saturated cells (i.e., S < 1), the effective edges lengths are computed using an upwind procedure: We further assume that ∆z n Figure 1. Schematic of the used staggered mesh.

Numerical Method
In order to avoid the development of a fully nonlinear system, the convective and the viscous terms in Equation (4) are discretized explicitly, while the velocity field and the hydraulic head in Equations (2) and (4) are discretized implicitly. A finite volume approximation of the continuity Equation (2) reads: Consistently, a finite difference approximation for the momentum Equations (4) reads: where Fu n i+ 1 2 ,k , Fw n i,k+ 1 2 contain all the explicit contributions of the convective and viscous terms, that in this paper are expressed through the explicit conservative formulation proposed in [39]. Substituting the discrete velocities u n+1 i± 1 2 ,k , w n+1 i,k± 1 2 as expressed in Equations (10) and (11) into Equation (9), we obtain a linear system for the unknown hydraulic head: where the term b n i,k contains all the known terms evaluated at the time t n : Equation (12) is a five-points diagonal system, symmetric and at least semi-positive definite (e.g., [36]). Therefore, it can be solved using a matrix-free implementation of the conjugate gradient method. Once the hydraulic head η n+1 is known, the velocity field can be easily updated through the explicit formulae in Equations (10) and (11).
As for the scalar transport, similarly to [38], we refer to a semi-implicit finite volume approximation based on upwind fluxes: where c n,up · is the upwind contribution defined as: If there is no sedimentation/erosion, then c n+1 i,k = c * i,k , otherwise the mass transfer is performed assuming that the fluid is in a local-static equilibrium. First, the dynamic part of the suspended sediment is computed through Equation (14), then the result is used to update explicitly the new sediment level through the discrete version of Equation (6): The sediment level at the time step n + 1 as resulting from Equation (17) allows for updating the volumes/edges lengths specified in Equation (8). Finally, the concentration at the time step n + 1 is updated according to the discrete version of the mass conservation law (7) that reads: We note that by substituting Equation (14) into (18) and integrating c n+1 i,k over the water column, one obtains the discrete mass conservation law for the suspended transport, where the source term is the discrete version of ∂z b ∂t (1 − φ), with z b total sediment level. However, the form in (14)- (18) is more general since it does not require a single value function z b = z b (x, t), but it can possibly be applied to cases with multiple sediment/water interfaces. For the chosen explicit discretization of the nonlinear convective terms, the method is stable under a Courant-Friedrichs-Lewy (CFL) type restriction based on the fluid velocity (e.g., [40]), The method becomes unconditionally stable if an Eulerian-Lagrangian scheme is adopted, in combination with a sub-step approach for the evolution of the concentration and sediment level [37,41].

Crank-Nicholson Time Discretization
For non-stationary problem we can improve the temporal accuracy by means of the so called θ-method. In order to do that, u n+1 and w n+1 in Equation (9) need to be substituted by u n+θ and w n+θ , while η n+1 in (10), (11) and (14) needs to be substituted by η n+θ . (u, v, η) n+θ are defined as: where θ is an implicit factor to be taken in the interval (0.5, 1] (e.g., [42]). The final system for the hydraulic head reads: where b n i,k becomes: We note that these modifications do not affect the structure of the linear system for the hydraulic head, since they are just rescaling it through a factor θ 2 .

Validation Tests
The numerical model was first validated against some standard numerical tests, for which reference solutions are available. In the first validation test, we considered the Blasius analytical solution for the laminar velocity profile in the boundary layer above a semi-infinite plate [43]. This validation test was obtained by simulating a plate 0.2 m long, assuming a uniform upstream velocity parallel to it and equal to u ∞ = 0.3 m/s. The two-dimensional x-z domain Ω = [−0.01, 0.2] × [0.0, 0.031] m was discretized according to a 600 × 300 grid, for a total of 180 000 elements having horizontal and vertical dimension of ∆x = 3.5 × 10 −4 m and ∆z = 1.03 × 10 −4 m, respectively. We set θ = 0.51, g = 1 m/s 2 , and ν = 10 −6 m 2 /s. As for the boundary conditions (BCs), we assumed the velocity u ∞ as the left BC, transmissive BCs at the right and at the top edges, and no-slip BC at the bottom plate, beginning from x = 0 in order to trigger the boundary layer.
The second validation test was performed considering the so called lid-driven cavity problem, which is another classical benchmark test [44]. The problem consists in a cavity Ω = [−0.5, 0.5] 2 m where the initial velocity field is (u, w) = 0.0 m. We set θ = 1, and g = 1 m/s 2 . We imposed (u, w) = (1, 0) m/s as the top BC and no-slip BCs at the other three boundaries. We repeated the test for two different values of the Reynolds number, Re = 400 and Re = 1 000, assuming the sediment level S ≡ 0. For both tests we discretized the domain according to a square grid 400 × 400, for a total of 160 000 elements having horizontal and vertical dimensions of ∆x = 2.5 × 10 −3 m and ∆z = 2.5 × 10 −3 m, respectively. This benchmark test was chosen due to its analogy to the inter-grain regions, where circulation cells develop bounded laterally and at the bottom by the grains and by the fine sediment, respectively.
The same lid-driven cavity test was performed introducing an erodible sediment at the bottom of the cavity, characterized by density ρ s = 1 553 kg/m 3 and porosity φ = 0.46. This resulted in considering mobile boundaries of the fluid domain due to the variation in time of the bed level z b , according to the following equation: where E and D are erosion and deposition rates, respectively. In Equation (23) we assumed D = 0, which corresponds to simulating a non-equilibrium erosive process. We considered the same domain Ω, grid discretization, BCs, and g as in the original lid-driven cavity test, but we filled the cavity with sediment up to 1/3 of the cavity height, set θ = 0.51, and performed the test for just one value of the Reynolds number, i.e., Re = 1 000. We defined the lowering of the fine sediment bed through the van Rijn erosion rate formula [45], which, expressed in [m/s], reads as follows: where ∆ is the fine sediment relative density , and T is the dimensionless excess of shear stress In the above formulas, ρ s is the fine sediment density, Θ is the Shields parameter, and Θ cr its critical value. The Shields parameter Θ is defined as: where τ b is the shear stress at the bottom, here defined as: To fasten the simulation, the erosion rate E from Equation (24) was multiplied by a factor 100. This test was introduced to validate the water and sediment mass conservation properties of the model, and to verify its robustness when dealing with time-varying changes in the boundaries of the fluid domain.

Numerical Experiments
Gravel bed rivers are often represented in laboratory experiments by covering the flume bed with spheres or hemispheres (e.g., [26,28,46,47]). We therefore assumed two simplified topographic configurations with the gravel represented by homogeneous spheres with different spacing, according to the most common simplified configurations used in laboratory flumes: in-line arrangement ( Figure 2a) and closest-packing arrangement (Figure 2b). In particular, in the first case, we examined the section where two consecutive spheres touch each other (Figure 2d), while in the second case we examined the section with the maximum inter-sphere spacing (Figure 2e).
To investigate the differences between a simplified and a real bed topography, we considered also a real gravel topography obtained from point-laser-scanning a longitudinal section in a laboratory flume covered by gravel (Figures 2c). In this case, we selected a 20 cm long section (Figure 2f), whose granulometric distribution was characterized by D 50 = 26.5 mm and D 90 = 29.5 mm. In order to set up equivalent geometry configurations, the two simplified cases described above were drawn considering spheres of 30.0 mm diameter.
In addition to the bed topography of the gravel matrix, in the numerical experiments we considered also the presence of inter-grain inerodible fine sediments at fixed level, to see the variations in the flow field depending on the depth. The fine sediment interface was schematized as a horizontal line, whose level was assigned according to four different filling rates: Z = 0, 0.25, 0.50, and 0.75, where Z = D/2 + z b D/2 and D is either the diameter of the spheres in the simplified case, or D 90 in the real topography. The vertical coordinate z was defined positive upwards, with z = 0 at the gravel crest. All the simulations were performed considering a 0.2 m long domain with a rigid top boundary at 0.03 m above the gravel crest. The domain was discretized according to a 400 × 600 grid for a total of 240 000 elements having horizontal and vertical dimension of dx = 5 × 10 −4 m and dz = 10 −4 m, respectively. We set θ = 1, and g = 9.81 m/s 2 , and ν = 10 −6 m 2 /s. The upstream and downstream boundary conditions were as follows: a uniform horizontal inflow velocity equal to 0.3 m/s and transmissive conditions downstream, with an assigned upstream/downstream pressure gradient equal to 2.5 × 10 −3 m/m. No slip boundary conditions were assigned at the bottom, and free slip conditions were assigned at the top boundary. We notice that solving the incompressible Navier-Stokes equations has the undeniable advantage of reducing the computational burden, while the simulation setup is equivalent to considering a free-surface flow with a bed slope equal to the pressure gradient.

Validation Tests
The Blasius solution is characterized by the following non linear third-order Ordinary Differential Equation (ODE): where f = u/u ∞ and ξ = z u ∞ 2νx is the so called Blasius coordinate. f is the primitive function of f , f is the second order derivative, f is the third order derivative. The reference solution is then obtained using a tenth-order Discontinuous Galerkin (DG) ODE solver [48]. In Figure 3a, the simulated velocity field is presented at time t = 10 s, and clearly shows the formation of the Blasius boundary layer. The horizontal velocity profile in the plane (x, ξ) and the velocity profile taken from the 25%, 50%, and 75% of the domain are compared with the exact Blasius solution in Figure 3b. A good agreement between the analytical solution and the numerical results can be observed, despite the small value of ν (i.e., 10 −6 m 2 /s). Furthermore, since the Blasius solution depends only on ξ, we expect the velocity profile to remain constant in the (x, ξ) plane, as is confirmed in Figure 3c.
As for the classical (i.e., fix bed) lid-driven cavity test, we calculated the velocity field in the cavity and compared it to the available reference solution given by [44]. Figures 4a and 5a show the velocity streamlines in the cavity at the final time t end = 100 s, while the comparison with the reference solution is shown in Figures 4b and 5b. In particular, these latter plots show the normal velocities passing through the lines {x = 0} and {z = 0}. For both the values of the Reynolds number, the figures show good agreement between the analytical solution and the numerical results, as it is also confirmed by the correct formation of secondary circulation cells at the two lower corners of the domain, as found also by [44].    The results of the same lid-driven cavity test, but performed on an erodible bed, are summarized in Figure 6. The figure shows the velocity field and the suspended sediment concentration field at times t = 20, 50, 100, 300 s, from which it is possible to appreciate that the erosion process is controlled by the main circulation cell, while the secondary circulation cells change in time according to the evolution of the water-sediment interface. The change in the suspended sediment concentration is due to the increase in the amount of sediments entrained from the bottom, as reflected by the progressive lowering of the bottom level in panels e-h. Since the cavity is a closed system, the progressive erosion of the bottom increases the total amount of suspended sediments in the domain.
In Figure 7 we report the measured mass exchange. The left side shows the time evolution of the total mass of suspended sediment and the total mass of deposited sediment, whereas the right side shows the balance between the free water mass and the mass of water constrained in the deposited sediment voids. Both sediment and water mass conservation are satisfied during the entire simulation with a precision determined by the tolerance used in the conjugate gradient algorithm (10 −10 in all cases shown here), hence possibly going down to the machine epsilon. Overall, this last validation test indicates the robustness of the model when dealing with sediment erosion and transport processes and with the associated time-varying evolution of the fluid domain boundaries. Model robustness is also indicated by the results obtained considering a 100 × 100 computational grid (10,000 elements having horizontal and vertical dimension of ∆x = 10 −2 m and ∆z = 10 −2 m; see Figure S1 in the Supplementary Materials), which are coherent with those shown in Figure 6, although in this latter case we clearly see a higher level of detail and lower numerical diffusion.

Numerical Experiments
In this section, we present the numerical results obtained for the simplified and real gravel bed topographies, that were simulated as representative cases of typical laboratory flume configurations. The simulated flow fields are shown in Figures 8-10 for all fine sediment filling rates-topography combinations. Results are presented in terms of horizontal velocity component and streamlines (left panels), and vorticity ω (right panels), where in a two-dimensional flow ω is defined as follows: y being the unit vector along the third dimension y. The fine grid resolution used to discretize the computational domain allowed to well capture the inter-grain flow structures, including the development of secondary circulations and stagnation points. In all cases, the flow field shown in Figures 8-10 refers to the end of the simulation (15 s), when the influence of the initial conditions was substantially lost and the flow field reached statistical steady state conditions with the development of periodic eddies generated by the variable topography ( Figure 11).  Despite advancements in experimental technologies exploited in recent studies (e.g., [22,24,25,32]), a detailed evaluation of the inter-grain flow field is still hampered by the small spatial scales of the geometries and processes at hand. However, being able to simulate the flow field at the inter-grain scale is certainly a desired goal, in that it allows for deriving the total shear stress distribution, which is the key quantity controlling fine sediment dynamics. According to the double-averaging approach, in gravel bed topographies, the total shear stress is expressed as the combination of viscous, turbulent, form-induced stresses, and form drag [23]. Quantifying the relative influence of the effective components (i.e., turbulent and form-induced stresses) and dispersive component (i.e., the form drag) is crucial to accurately estimate the entrainment and transport of the inter-grain fine sediment. In this regard, the above results suggest that the proposed numerical model has a great potential for applications in the context of fine sediment transport dynamics in gravel bed rivers. Among the major advantages is the possibility to acquire the fine resolution flow field in continuous, whereas experimental measurements are typically performed at discrete positions that get rarer going further in depth below the gravel crest. This means loosing information that a continuous numerical measurement can ensure, such as the existence of the inter-grain secondary circulations observed in Figures 8-10. Such structures determine a double inflection in the velocity profile in all topographic configurations for decreasing fine sediment filling rates. This clearly emerges from the vertical profiles of the horizontal component of the velocity shown in Figure 12, at chosen sections representative of the gravel crest and of the inter-grain cavity (see Figure 2). The presented profiles are averaged along 5 s of the simulation (from 10 to 15 s, with an output resolution of 0.1 s), for the same fine sediment filling rates-topography combinations of Figures 8-10. The instantaneous velocity profiles used to compute the averaged profiles are also shown for the inter-grain cavity section (thin continuous lines), which indicate the high non-stationary behavior of the flow field in the roughness layer [23]. The same figure but for the vertical component of the velocity is shown in Figure 13. As expected, both Figures 12 and 13 show clear differences depending on fine sediment filling rates and gravel bed topography. The bed topography influence on the flow field confirms that special attention should be paid when using simplified gravel bed geometries in laboratory flume experiments. In fact, despite they offer undeniable advantages in terms of simplification of the experimental setup (e.g., by allowing for an easier definition of porosity, granulometric distribution, and bed topography), a spheres-covered bed may not be entirely representative of a natural water-worked gravel bed due to different macro-roughness structures. Strong differences are particularly evident for the vertical component of the velocity (Figure 13), which is responsible for lifting (i.e., eroding) fine material at the bottom [22,24,25]. In this regard, we note that the instantaneous profiles of vertical velocity undergo large changes from negative (i.e., downward velocity) to positive (i.e., upward velocity) values in all topographic configurations. Time variability in the velocity field is also evident from the profiles of the horizontal component (Figure 12), indicating the existence of highly non-stationary inter-grain recirculation cells. This is confirmed also in the vorticity field and particle tracking videos available in the Supplementary Materials.

Conclusions
A second order semi-implicit numerical scheme on staggered Cartesian meshes for the incompressible Navier-Stokes equations, based on the method proposed by [36][37][38] in presence of a time dependent sedimentation/erosion process, was derived. In the scheme, we defined the hydraulic head in the cells centers and the velocity at the cells interfaces. By formally substituting the discrete momentum equations into the discrete continuity equation, we obtained a symmetric semi-positive definite linear system where the only unknown is the hydraulic head at the new time step. The system is then solved using a fast iterative linear solver such as the conjugate gradient algorithm [49]. We note that the method is built in such a way that the computation in each cell involves only its direct neighbors. This makes the algorithm particularly suitable to parallelization, since the data that need to be synchronized are limited to the single layer of cells surrounding each parallel region.
For the entrainment and deposition of the sediment we used an explicit finite volume scheme in combination with a general mass flux between suspended and deposited sediment. The deposition of the suspended sediment changes the effective domain sizes in terms of volume and edges length in the cells. The method is mass-conservative and limited in the time discretization by a classical CFL-type time restriction based on the local fluid velocity. However, if the convective-viscous terms are solved by an Eulerian-Lagrangian method combined with a local time stepping/subcycling approach for the sediment dynamics, the method becomes unconditionally stable. Furthermore, compared to [36], the pressurized system allows for avoiding the solution of the mildly nonlinear contribution through the Nested-Newton approach [50], which ultimately reduces the computational time thus allowing for a fine resolution in the mesh.
The method was validated against some classical benchmarks, i.e., the Blasius boundary layer and the lid-driven cavity test. Moreover, a modified version of the lid-driven cavity test was run, to verify the conservation of sediment and water mass in presence of erodible sediment, and the robustness of the model in presence of time-varying boundaries of the fluid domain.
Once validated, the model was used to simulate three simple cases representative of typical experiments for gravel bed rivers at different filling rates. The numerical results for the inter-grain flow field show the formation of main and secondary circulation cells, forced by the presence of the macro-roughness elements, which generates a double inflection in the time-averaged velocity profile for the lower filling rates. This information would probably be lost if performing only experimental measurements of the flow dynamics, which are possible just at discrete points, especially below the gravel crest. Furthermore, the use of a numerical model simplifies the repetition of the experiments considering different topographies, which contributes at improving the understanding of the geometry role in the stresses distribution.
We note that the model is two-dimensional, therefore it does not account for the three-dimensional effects and its results are not immediately representative of the real case. While a full description of the inter-grain flow field and turbulence structure would require the use of a three-dimensional model coupled with a proper turbulence closure scheme, even in this form the proposed model provides useful clues on the approximation effects introduced when using simplified geometries to represent real topographies.
Future work will concern the extension of the present approach to the complete VOF (Volume Of Fluid) method such as proposed in [36], and its extension to the three-dimensions in the presence of erodible sediment, together with the inclusion of a proper turbulence closure and high-performance parallelization standards. Funding: Part of this work was supported by the projects Sediplan-r (FESR1002) and LTFD Laboratory of Thermo Fluid Dynamics (FESR1029), financed by the European Regional Development Fund (ERDF) Investment for Growth and Jobs Programme 2014-2020, and CRC project HM: Hydropeaking mitigation, financed by the Free University of Bozen-Bolzano.