*3.1. Validation Tests*

The Blasius solution is characterized by the following non linear third-order Ordinary Differential Equation (ODE):

$$f'''' + f \cdot f'' = 0 \qquad f(0) = 0, \quad f'(0) = 0, \quad \lim\_{\substack{\mathbb{F} \to \infty}} f'(\xi) = 1,\tag{27}$$

where *f* = *u*/*u*<sup>∞</sup> and *ξ* = *z u*<sup>∞</sup> <sup>2</sup>*ν<sup>x</sup>* is the so called Blasius coordinate. *<sup>f</sup>* 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.

**Figure 3.** (**a**) Field map of the horizontal velocity component *u*, showing the formation of the Blasius boundary layer; (**b**) comparison between the numerical solution and the analytical Blasius solution at three different locations of the domain; and (**c**) magnitude of the horizontal velocity component *u* in the (*x*, *ξ*) plane.

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−<sup>6</sup> m2/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 *tend* = 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].

**Figure 4.** (**a**) Streamlines in the cavity at the final time *t* = 100 s; and (**b**) comparison with the reference solution [44] for *Re* = 400.

**Figure 5.** (**a**) Streamlines in the cavity at the final time *t* = 100 s; and (**b**) comparison with the reference solution [44] for *Re* = 1 000.

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−<sup>10</sup> 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−<sup>2</sup> m and Δ*z* = 10−<sup>2</sup> 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.

**Figure 6.** (**a**–**d**) Evolution of the streamlines in the cavity with the erodible bed at times *t* = 20, 50, 100, 300 s and coloured velocity magnitude |*v*|, where *v* denotes the velocity vector; and (**e**–**h**) volume concentration of suspended sediment.

**Figure 7.** Time evolution of (**a**) sediment and (**b**) water mass, rescaled with respect to the initial total mass.

#### *3.2. 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:

$$
\vec{\omega} = \left(\frac{\partial w}{\partial \mathbf{x}} - \frac{\partial u}{\partial z}\right) \vec{y} \,\tag{28}
$$

*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).

**Figure 8.** Horizontal velocity component *u* and streamlines (**left panels**), and magnitude of the vorticity |*ω*-| (**right panels**) for the spheres in-line arrangement.

**Figure 9.** Horizontal velocity component *u* and streamlines (**left panels**), and magnitude of the vorticity |*ω*-| (**right panels**) for the spheres closest-packing arrangement.

*Water* **2020**, *12*, 690

**Figure 10.** Horizontal velocity component *u* and streamlines (**left panels**), and magnitude of the vorticity |*ω*-| (**right panels**) for the real gravel topography configuration.

**Figure 11.** Time evolution of the horizontal (*u*) and vertical (*w*) velocity component intensities, averaged over the region below the gravel crest level, and relative to the real gravel bed configuration with fine sediment filling rate *Z* = 0.5.

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.

**Figure 12.** Vertical profiles of the horizontal velocity component *u* at chosen sections representative of the gravel crest and of the inter-grain cavity (see Figure 2), for the in-line arrangement (**left panels**), closest-packing arrangement (**central panels**), and real gravel bed configuration (**right panels**). Thick lines indicate profiles averaged from 10 to 15 s of simulation, while thin lines indicate instantaneous profiles every 0.1 s.

**Figure 13.** Vertical profiles of the vertical velocity component *w* at chosen sections representative of the gravel crest and of the inter-grain cavity (see Figure 2), for the in-line arrangement (**left panels**), closest-packing arrangement (**central panels**), and real gravel bed configuration (**right panels**). Thick lines indicate profiles averaged from 10 to 15 s of simulation, while thin lines indicate instantaneous profiles every 0.1 s.

#### **4. 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–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.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/3/690/s1, Figure S1: the numerical results for the lid-driven cavity test with erodible bed obtained using a coarser domain discretization; video S1: the vorticity field for the closest-packing arrangement; video S2: particle tracking for the closest-packing arrangement.

**Author Contributions:** Conceptualization, M.T., S.P. and G.S.; methodology, M.T., S.P., and G.S.; validation, M.T.; formal analysis, M.T., S.P. and G.S.; data curation M.T. and S.P.; writing–original draft preparation, G.S. and S.P.; writing–review and editing, all authors; supervision, M.R. All authors have read and agreed to the published version of the manuscript.

**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.

**Acknowledgments:** We are grateful to E. Spilone for relevant discussion during the preparation of this work.

**Conflicts of Interest:** The authors declare no conflict of interest.
