Next Article in Journal
Detection of Helminth Ova in Wastewater Using Recombinase Polymerase Amplification Coupled to Lateral Flow Strips
Next Article in Special Issue
Seasonal Precipitation Variability and Gully Erosion in Southeastern USA
Previous Article in Journal
Groundwater Contribution to Sewer Network Baseflow in an Urban Catchment-Case Study of Pin Sec Catchment, Nantes, France
Previous Article in Special Issue
A Fuzzy Transformation of the Classic Stream Sediment Transport Formula of Yang
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Faculty of Science and Technology, Free University of Bozen-Bolzano, Universitätsplatz - Piazza Università 5, 39100 Bozen-Bolzano, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Water 2020, 12(3), 690; https://doi.org/10.3390/w12030690
Submission received: 14 December 2019 / Revised: 15 February 2020 / Accepted: 18 February 2020 / Published: 3 March 2020
(This article belongs to the Special Issue Modeling of Soil Erosion and Sediment Transport)

Abstract

:
The selective trapping and erosion of fine 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 fine sediment erosion and transport in gravel bed rivers, a comprehensive overview of the leading processes is hampered by the difficulties in performing flow field 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 efficient, and allows for a fine 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.

1. 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].
Occasionally, armored stream beds may be subject to high sediment loads, for example during dam flushing or removal operations, or during flood events associated to large sediment input from lateral inflows. Under these conditions, finer particles infiltrate into the void spaces of the immobile coarse bed grains, according to a selective trapping mechanism that is reciprocal to the selective erosion that caused bed armoring [11]. If the infiltration of fine particles into the coarse bed interstices is extensive, the volume of voids drops increasing the compactness of the stream bed texture, thus decreasing its hydraulic conductivity and increasing its resistance to flow (e.g., [12,13]). This process is known as colmation or clogging (see, for example, [14] for a review), and has significant impacts on stream ecology (e.g., [15,16,17,18]), exchanges of water, dissolved substances and heat with the underlying hyporheic zone and groundwater [19], and flow and turbulence structure [20,21,22]. Under high flow conditions, the armour layer can break up and the entire river bed becomes mobilized, hence resetting the bed morphology and grain size distribution. However, if the erosion capacity of the flow is not sufficient to remove the coarse grains, as soon as the sediment load from upstream declines, the reestablishment of sediment supply-limited conditions reactivates selective erosion of finer particles, sustaining the armoring of the stream bed.
The cleaning dynamics controlling the erosion of finer particles from coarse granular beds is inherently different from that typical of uniform sized beds. In fact, in armored beds, the presence of macro-roughness due to the coarser particles alters the flow structure and, consequently, the distribution of the stress components below the gravel crest level. In these beds, besides turbulent stresses, form-induced stresses and form drag also contribute to the total shear stress distribution [23]. In addition, the vertical component of the stress responsible for lifting and transporting fine material was found to be decreasing below gravel crest level [22,24,25]. This alteration is mirrored in the reduced sediment entrainment and transport capacity of the flow, which is affected also by the smaller fine sediment-water active interface with respect to that of a uniform bed.
It is therefore clear that the traditional formulae derived for uniform bed cases fail to describe erosion and sediment transport processes over immobile gravel beds. In fact, these formulae do not account for the reduction in the effective part of the shear stress, nor the reduction in the fine sediment-water active interface. In this respect, performing laboratory experiments is a common methodological practice for investigating selective transport dynamics in armored beds and gaining useful elements to derive empirical formulas of fine sediment transport (e.g., [26,27,28,29]). Experimental research is typically carried out in laboratory flumes using laser-scanner (e.g., [29]) or digital photogrammetry (e.g., [30]) to measure the changes in the topography. The track of these changes in the fine sediment level inside the gravel matrix is usually coupled with measures of transport rate (trough sieves or density cells (e.g., [27,29]) or concentration (e.g., [26,31]) to quantify the fine sediment transport and/or erosion rate between the gravel. Based on these experimental approaches, useful fine sediment transport formulae has been proposed in previous literature, as for instance in [26,27,28,29].
Despite the above cited empirical formulae and despite the examples of direct measurements of the flow field in the roughness layer [23] (e.g., [22,24,25,32]), a comprehensive framework on the inter-grain flow and sediments dynamics in gravel bed rivers still poses some scientific challenges. These challenges are primarily due to the operational difficulties to perform velocity measurements far below the gravel crest level and to quantify the relative contribution of the form drag to the total shear stress [25]. In addition, a fair comparison among existing studies is not obvious due to the differences in the bed topographies, which chiefly controls the distribution of the shear stress components [25], thus hampering the derivation of general considerations. In this context, fine-scale numerical models can offer a valid alternative to overcome the inherent difficulties of fine-resolution, inter-grain experimental measurements. At the same time, they can easy the investigation of the role of the geometry in affecting the stresses distribution, provided that the setup and repetition of laboratory experiments with different configurations is not a minor matter. While examples of Direct Numerical 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.

2. Materials and Methods

2.1. Numerical Model

2.1.1. Governing Equations

The governing equations are the two-dimensional incompressible Navier-Stokes, that are composed by the momentum conservation:
u t + u u x + u w z = p x + ν 2 u x 2 + 2 u z 2 , w t + w u x + w w z = p z + ν 2 w x 2 + 2 w z 2 g ,
and the incompressibility condition:
u x + w z = 0 ,
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 + u u x + u w z = g η x + ν 2 u x 2 + 2 u z 2 , w t + w u x + w w z = g η z + ν 2 w x 2 + 2 w z 2 .
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:
c t + u c x + ( w w c ) c z = f c S ( S , c )
S t = f S c ( S , c ) ,
where w c is the falling velocity of the suspended sediment; f S c and f c S 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 S c and f c S . 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:
d d t c + S ( 1 ϕ ) = 0 .

2.1.2. 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 Ω i , k , where Ω i , k = [ x i 1 2 , x i + 1 2 ] × [ z k 1 2 , z k + 1 2 ] are rectangular control volumes centered in ( x i , z k ) . The mesh is assumed to be homogeneous with Δ x = x i + 1 2 x i 1 2 and Δ z = z k + 1 2 z k 1 2 , for i = 1 , , I m a x and k = 1 , , K m a x . 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 η = η i k , c = c i k , S = S i k . 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 i k ) . 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 i ± 1 2 , k n and Δ x i , k ± 1 2 n , see Figure 1. For non-saturated cells (i.e., S < 1 ), the effective edges lengths are computed using an upwind procedure:
Δ z i + 1 2 , k n = V i , k n Δ x u i + 1 2 , k n > 0 V i + 1 , k n Δ x u i + 1 2 , k n < 0 max ( V i , k n , V i + 1 , k n ) Δ x u i + 1 2 , k n = 0 Δ x i , k + 1 2 n = 0 S i , k = 1 Δ x S i , k 1
We further assume that Δ z i ± 1 2 , k n = Δ x i , k ± 1 2 n = 0 for any ( i , k ) where S i , k = 1 .

2.1.3. 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:
Δ z i + 1 2 , k n u i + 1 2 , k n + 1 Δ z i 1 2 , k n u i 1 2 , k n + 1 + Δ x i , k + 1 2 n w i , k + 1 2 n + 1 Δ x i , k 1 2 n w i , k 1 2 n + 1 = 0 .
Consistently, a finite difference approximation for the momentum Equations (4) reads:
u i + 1 2 , k n + 1 = F u i + 1 2 , k n Δ t g η i + 1 , k n + 1 η i , k n + 1 Δ x
w i , k + 1 2 n + 1 = F w i , k + 1 2 n Δ t g η i , k + 1 n + 1 η i , k n + 1 Δ z ,
where F u i + 1 2 , k n , F w i , k + 1 2 n 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 i ± 1 2 , k n + 1 , w i , k ± 1 2 n + 1 as expressed in Equations (10) and (11) into Equation (9), we obtain a linear system for the unknown hydraulic head:
Δ t g Δ x Δ z i + 1 2 , k n η i + 1 , k n + 1 η i , k n + 1 Δ z i 1 2 , k n η i , k n + 1 η i 1 , k n + 1 + Δ t g Δ z Δ x i , k + 1 2 n η i , k + 1 n + 1 η i , k n + 1 Δ x i , k 1 2 n η i , k n + 1 η i , k 1 n + 1 = b i , k n ,
where the term b i , k n contains all the known terms evaluated at the time t n :
b i , k n = Δ z i + 1 2 , k n F u i + 1 2 , k n Δ z i 1 2 , k n F u i 1 2 , k n + Δ x i , k + 1 2 n F w i , k + 1 2 n Δ x i , k 1 2 n F w i , k 1 2 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:
c i , k * = c n + Δ t Δ x Δ z Δ z i + 1 2 , k n c i + 1 2 , k n , u p u i + 1 2 , k n + 1 Δ z i 1 2 , k n c i 1 2 , k n , u p u i 1 2 , k n + 1 Δ x i , k + 1 2 n c i , k + 1 2 n , u p w i , k + 1 2 n + 1 w s Δ x i , k 1 2 n c i , k 1 2 n , u p w i , k 1 2 n + 1 w s ,
where c · n , u p is the upwind contribution defined as:
c i + 1 2 , k n , u p ( V ) = 1 2 c i , k n ( | V | + V ) + c i + 1 , k n ( | V | V )
c i , k + 1 2 n , u p ( V ) = 1 2 c i , k n ( | V | + V ) + c i , k + 1 n ( | V | V )
If there is no sedimentation/erosion, then c i , k n + 1 = 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):
S i , k n + 1 = S i , k n + Δ t Δ x Δ z Q S C ( S i , k n , c i , k * ) , where Q S C ( S i , k n , c i , k * ) = Ω i , k f S C ( S i , k n , c i , k * ) .
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:
c i , k n + 1 = c i , k * + S i , k n S i , k n + 1 ( 1 ϕ ) .
We note that by substituting Equation (14) into (18) and integrating c i , k n + 1 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]),
Δ t 1 | u | Δ x + | w | Δ z + 2 ν Δ x 2 + 2 ν Δ z 2 .
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].

2.1.4. 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:
( u , v , η ) n + θ = θ · ( u , v , η ) n + 1 + ( 1 θ ) · ( u , v , η ) n ,
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:
Δ t g θ 2 Δ x Δ z i + 1 2 , k n η i + 1 , k n + 1 η i , k n + 1 Δ z i 1 2 , k n η i , k n + 1 η i 1 , k n + 1 + Δ t g θ 2 Δ z Δ x i , k + 1 2 n η i , k + 1 n + 1 η i , k n + 1 Δ x i , k 1 2 n η i , k n + 1 η i , k 1 n + 1 = b i , k n
where b i , k n becomes:
b i , k n = θ Δ z i + 1 2 , k n F u i + 1 2 , k n Δ z i 1 2 , k n F u i 1 2 , k n + Δ x i , k + 1 2 n F w i , k + 1 2 n Δ z i 1 2 , k n F w i , k 1 2 n + θ ( 1 θ ) Δ t g Δ x Δ z i + 1 2 , k n η i + 1 , k n η i , k n Δ z i 1 2 , k n η i , k n η i 1 , k n + θ ( 1 θ ) Δ t g Δ z Δ x i , k + 1 2 n η i , k + 1 n η i , k n Δ x i , k 1 2 n η i , k n η i , k 1 n ( 1 θ ) Δ z i + 1 2 , k n u i + 1 2 , k n Δ z i 1 2 , k n u i 1 2 , k n + Δ x i , k + 1 2 n w i , k + 1 2 n Δ z i , k 1 2 n w i , k 1 2 n .
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 .

2.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, R e = 400 and R e = 1000 , 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 = 1553 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:
z b t = D E ( 1 ϕ ) ,
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., R e = 1000 . 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:
E = 0.00033 g Δ d s d * 0.3 T 1.5 ,
where Δ is the fine sediment relative density Δ = ρ s ρ 1 , d s is the fine sediment median diameter, d * is the dimensionless grain size d * = d s Δ g ν 2 1 / 3 , and T is the dimensionless excess of shear stress T = Θ Θ c r Θ c r . In the above formulas, ρ s is the fine sediment density, Θ is the Shields parameter, and Θ c r its critical value. The Shields parameter Θ is defined as:
Θ = τ b Δ ρ g d s ,
where τ b is the shear stress at the bottom, here defined as:
τ b = μ u ( z ) z z = z b .
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 (Figure 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 d x = 5 × 10 4 m and d z = 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.

3. Results and Discussion

3.1. Validation Tests

The Blasius solution is characterized by the following non linear third-order Ordinary Differential Equation (ODE):
f + f · f = 0 f ( 0 ) = 0 , f ( 0 ) = 0 , lim ξ f ( ξ ) = 1 ,
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]. Figure 4a and Figure 5a show the velocity streamlines in the cavity at the final time t e n d = 100 s, while the comparison with the reference solution is shown in Figure 4b and Figure 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.

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 Figure 8, Figure 9 and Figure 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:
ω = w x u z y ,
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 Figure 8, Figure 9 and Figure 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 Figure 8, Figure 9 and Figure 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 Figure 8, Figure 9 and Figure 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 Figure 12 and Figure 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.

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

Supplementary Materials

The following are available online at https://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.

References

  1. Schleiss, A.J.; Franca, M.J.; Juez, C.; Cesare, G.D. Reservoir sedimentation. J. Hydraul. Res. 2016, 54, 595–614. [Google Scholar] [CrossRef]
  2. Kaffas, K.; Hrissanthou, V.; Sevastas, S. Modeling hydromorphological processes in a mountainous basin using a composite mathematical model and ArcSWAT. CATENA 2018, 162, 108–129. [Google Scholar] [CrossRef]
  3. ICOLD. Sedimentation and Sustainable Use of Reservoirs and River Systems; Technical Report; International Commission on Large Dams: Paris, France, 2009. [Google Scholar]
  4. Williams, G.P.; Wolman Gordon, M. Downstream Effects of Dams on Alluvial Rivers; Geological Survey Professional Paper 1286; U.S. Government Printing Office: Washington, DC, USA, 1984. [CrossRef] [Green Version]
  5. Brandt, S. Classification of geomorphological effects downstream of dams. CATENA 2000, 40, 375–401. [Google Scholar] [CrossRef]
  6. Juez, C.; Hassan, M.A.; Franca, M.J. The Origin of Fine Sediment Determines the Observations of Suspended Sediment Fluxes Under Unsteady Flow Conditions. Water Resour. Res. 2018, 54, 5654–5669. [Google Scholar] [CrossRef]
  7. Kondolf, G.M. Hungry Water: Effects of Dams and Gravel Mining on River Channels. Environ. Manag. 1997, 21, 533–551. [Google Scholar] [CrossRef]
  8. Shen, H.W.; Lu, J. Development and Prediction of Bed Armoring. J. Hydraul. Eng. 1983, 109, 611–629. [Google Scholar] [CrossRef]
  9. Dietrich, W.; Kirchner, J.; Ikeda, H.; Iseya, F. Sediment Supply and Development of Coarse Surface Layer in Gravel Bedded Rivers. Nature 1989, 340. [Google Scholar] [CrossRef]
  10. Wilcock, P.R.; DeTemple, B.T. Persistence of armor layers in gravel-bed streams. Geophys. Res. Lett. 2005, 32. [Google Scholar] [CrossRef] [Green Version]
  11. Bui, V.; Bui, M.; Rutschmann, P. Advanced Numerical Modeling of Sediment Transport in Gravel-Bed Rivers. Water 2019, 11, 550. [Google Scholar] [CrossRef] [Green Version]
  12. Schälchli, U. The clogging of coarse gravel river beds by fine sediment. Hydrobiologia 1992, 235–236, 189–197. [Google Scholar] [CrossRef]
  13. Wu, F.C.; Huang, H.T. Hydraulic Resistance Induced by Deposition of Sediment in Porous Medium. J. Hydraul. Eng. 2000, 126. [Google Scholar] [CrossRef]
  14. Wharton, G.; Mohajeri, S.H.; Righetti, M. The pernicious problem of streambed colmation: A multi-disciplinary reflection on the mechanisms, causes, impacts, and management challenges. Wiley Interdiscip. Rev. Water 2017, 4, e1231. [Google Scholar] [CrossRef]
  15. Brunke, M.; Gonser, T. The ecological significance of exchange processes between rivers and groundwater. Freshw. Biol. 1997, 37, 1–33. [Google Scholar] [CrossRef] [Green Version]
  16. Kemp, P.; Sear, D.; Collins, A.; Naden, P.; Jones, I. The impacts of fine sediment on riverine fish. Hydrol. Process. 2011, 25, 1800–1821. [Google Scholar] [CrossRef]
  17. Jones, J.; Collins, A.; Naden, P.; Sear, D. The relationship between fine sediment and macrophytes in rivers. River Res. Appl. 2012, 28, 1006–1018. [Google Scholar] [CrossRef]
  18. Jones, J.I.; Murphy, J.F.; Collins, A.L.; Sear, D.A.; Naden, P.S.; Armitage, P.D. The impact of fine sediment on macro-invertebrates. River Res. Appl. 2012, 28, 1055–1071. [Google Scholar] [CrossRef]
  19. Krause, S.; Hannah, D.M.; Fleckenstein, J.H. Hyporheic hydrology: Interactions at the groundwater-surface water interface. Hydrol. Process. 2009, 23, 2103–2107. [Google Scholar] [CrossRef]
  20. Sambrook Smith, G.H.; Nicholas, A.P. Effect on flow structure of sand deposition on a gravel bed: Results from a two-dimensional flume experiment. Water Resour. Res. 2005, 41. [Google Scholar] [CrossRef] [Green Version]
  21. Wren, D.; Langendoen, E.; Kuhnle, R. Effects of sand addition on turbulent flow over an immobile gravel bed. J. Geophys. Res. Earth Surf. 2011, 116, 1–12. [Google Scholar] [CrossRef] [Green Version]
  22. Mohajeri, S.H.; Righetti, M.; Wharton, G.; Romano, G.P. On the structure of gravel-bed flow with intermediate submergence: Implications for sediment transport. Adv. Water Resour. 2016, 92, 90–104. [Google Scholar] [CrossRef]
  23. Nikora, V.; Goring, D.; McEwan, I.; Griffiths, G. Spatially averaged open-channel flow over rough bed. J. Hydraul. Eng. 2001, 127, 123–133. [Google Scholar] [CrossRef]
  24. Mignot, E.; Barthelemy, E.; Hurter, D. Double-averaging analysis and local flow chracterization of near-bed turbulence in gravel-bed channel flow. J. Fluid Mech. 2009, 618, 279–303. [Google Scholar] [CrossRef]
  25. Dey, S.; Das, R. Gravel-bed hydrodynamics: Double-averaging approach. J. Hydraul. Eng. 2012, 138, 707–725. [Google Scholar] [CrossRef]
  26. Grams, P.E.; Wilcock, P.R. Equilibrium entrainment of fine sediment over a coarse immobile bed. Water Resour. Res. 2007, 43. [Google Scholar] [CrossRef] [Green Version]
  27. Kuhnle, R.A.; Wren, D.G.; Langendoen, E.J.; Rigby, J.R. Sand transport over an immobile gravel bed substrate. J. Hydraul. Eng. 2013, 139, 167–176. [Google Scholar] [CrossRef] [Green Version]
  28. Grams, P.E.; Wilcock, P.R. Transport of fine sediment over a coarse, immobile riverbed. J. Geophys. Res. Earth Surf. 2015, 119, 188–211. [Google Scholar] [CrossRef]
  29. Kuhnle, R.A.; Langendoen, E.J.; Wren, D.G. Prediction of sand transport over immobile gravel from supply-limited to capacity conditions. J. Hydraul. Eng. 2017, 143. [Google Scholar] [CrossRef]
  30. Bertin, S.; Friedrich, H. Effects of Sand Addition and Bed Flushing on Gravel Bed Surface Microtopography and Roughness. Water Resour. Res. 2019. [Google Scholar] [CrossRef]
  31. Kuhnle, R.A.; Wren, D.G.; Langendoen, E.J. Erosion of Sand from a Gravel Bed. J. Hydraul. Eng. 2016, 142, 04015052. [Google Scholar] [CrossRef]
  32. Wren, D.G.; Kuhnle, R.A.; Langendoen, E.J.; Rigby, J.R. Turbulent Flow and Sand Transport over a Cobble Bed in a Laboratory Flume. J. Hydraul. Eng. 2014, 140, 04014001. [Google Scholar] [CrossRef]
  33. Fornarelli, F.; Vittori, G. Oscillatory boundary layer close to a rough wall. Eur. J. Mech. B/Fluids 2009, 28, 283–295. [Google Scholar] [CrossRef]
  34. Ghodke, C.D.; Apte, S.V. DNS study of particle-bed–turbulence interactions in an oscillatory wall-bounded flow. J. Fluid Mech. 2016, 792, 232–251. [Google Scholar] [CrossRef]
  35. Mazzuoli, M.; Blondeaux, P.; Simeonov, J.; Calantoni, J. Direct numerical simulation of the oscillatory flow around a sphere resting on a rough bottom. J. Fluid Mech. 2017, 822, 235–266. [Google Scholar] [CrossRef] [Green Version]
  36. Casulli, V. A semi-implicit numerical method for the free-surface Navier–Stokes equations. Int. J. Numer. Methods Fluids 2014, 74, 605–622. [Google Scholar] [CrossRef]
  37. Casulli, V.; Zanolli, P. Semi-implicit numerical modeling of nonhydrostatic free-surface flows for environmental problems. Math. Comput. Model. 2002, 36, 1131–1149. [Google Scholar] [CrossRef]
  38. Casulli, V.; Zanolli, P. High resolution methods for multidimensional advection–diffusion problems in free-surface hydrodynamics. Ocean Model. 2005, 10, 137–151. [Google Scholar] [CrossRef]
  39. Stelling, G.S.; Duinmeijer, S.P.A. A staggered conservative scheme for every Froude number in rapidly varied shallow water flows. Int. J. Numer. Methods Fluids 2003, 43, 1329–1354. [Google Scholar] [CrossRef]
  40. Gross, E.S.; Bonaventura, L.; Rosatti, G. Consistency with continuity in conservative advection schemes for free-surface models. Int. J. Numer. Methods Fluids 2002, 38, 307–327. [Google Scholar] [CrossRef]
  41. Casulli, V. Eulerian-Lagrangian methods for the Navier-Stokes equations at high Reynolds number. Int. J. Numer. Methods Fluids 1988, 8, 1349–1360. [Google Scholar] [CrossRef]
  42. Casulli, V.; Cattani, E. Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow. Comput. Math. Appl. 1994, 27, 99–112. [Google Scholar] [CrossRef] [Green Version]
  43. Blasius, H. Grenzschichten in Flüssigkeiten mit kleiner Reibung. Z. Math. Phys. 1908, 56, 1–37. [Google Scholar]
  44. Ghia, U.; Ghia, K.; Shin, C.T. High-resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J. Comput. Phys. 1982, 48, 387–411. [Google Scholar] [CrossRef]
  45. van Rijn, L.C. Sediment pick-up functions. J. Hydraul. Eng. 1984, 110, 1494–1502. [Google Scholar] [CrossRef]
  46. Manes, C.; Pokrajac, D.; McEwan, I.; Nikora, V. Turbulence structure of open channel flows over permeable and impermeable beds: A comparative study. Phys. Fluids 2009, 21, 125109. [Google Scholar] [CrossRef] [Green Version]
  47. Yager, E.M.; Kirchner, J.W.; Dietrich, W.E. Calculating bed load transport in steep boulder bed channels. Water Resour. Res. 2007, 43. [Google Scholar] [CrossRef]
  48. Dumbser, M. Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier–Stokes equations. Comput. Fluids 2010, 39, 60–76. [Google Scholar] [CrossRef]
  49. Hestenes, M.R.; Stiefel, E. Methods of conjugate gradients for solving linear systems. J. Res. Natl. Bur. Stand. 1952, 49, 409–436. [Google Scholar] [CrossRef]
  50. Casulli, V.; Zanolli, P. Iterative solutions of mildly nonlinear systems. J. Comput. Appl. Math. 2012, 236, 3937–3947. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic of the used staggered mesh.
Figure 1. Schematic of the used staggered mesh.
Water 12 00690 g001
Figure 2. Schematic of the simulated gravel bed topographies: (a) in-line arrangement, view from above; (b) closest-packing arrangement, view from above; (c) real gravel bed, 3D view; (d) in-line arrangement, longitudinal section; (e) closest-packing arrangement, longitudinal section; (f) real gravel bed, longitudinal section. In panels (df), the horizontal dashed lines indicate the fine sediment filling rate, while the vertical lines indicate the gravel crest (dotted) and inter-grain cavity (continuous) sections.
Figure 2. Schematic of the simulated gravel bed topographies: (a) in-line arrangement, view from above; (b) closest-packing arrangement, view from above; (c) real gravel bed, 3D view; (d) in-line arrangement, longitudinal section; (e) closest-packing arrangement, longitudinal section; (f) real gravel bed, longitudinal section. In panels (df), the horizontal dashed lines indicate the fine sediment filling rate, while the vertical lines indicate the gravel crest (dotted) and inter-grain cavity (continuous) sections.
Water 12 00690 g002
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.
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.
Water 12 00690 g003
Figure 4. (a) Streamlines in the cavity at the final time t = 100 s; and (b) comparison with the reference solution [44] for R e = 400 .
Figure 4. (a) Streamlines in the cavity at the final time t = 100 s; and (b) comparison with the reference solution [44] for R e = 400 .
Water 12 00690 g004
Figure 5. (a) Streamlines in the cavity at the final time t = 100 s; and (b) comparison with the reference solution [44] for R e = 1 000 .
Figure 5. (a) Streamlines in the cavity at the final time t = 100 s; and (b) comparison with the reference solution [44] for R e = 1 000 .
Water 12 00690 g005
Figure 6. (ad) 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 (eh) volume concentration of suspended sediment.
Figure 6. (ad) 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 (eh) volume concentration of suspended sediment.
Water 12 00690 g006
Figure 7. Time evolution of (a) sediment and (b) water mass, rescaled with respect to the initial total mass.
Figure 7. Time evolution of (a) sediment and (b) water mass, rescaled with respect to the initial total mass.
Water 12 00690 g007
Figure 8. Horizontal velocity component u and streamlines (left panels), and magnitude of the vorticity ω (right panels) for the spheres in-line arrangement.
Figure 8. Horizontal velocity component u and streamlines (left panels), and magnitude of the vorticity ω (right panels) for the spheres in-line arrangement.
Water 12 00690 g008
Figure 9. Horizontal velocity component u and streamlines (left panels), and magnitude of the vorticity ω (right panels) for the spheres closest-packing 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 12 00690 g009
Figure 10. Horizontal velocity component u and streamlines (left panels), and magnitude of the vorticity ω (right panels) for the real gravel topography configuration.
Figure 10. Horizontal velocity component u and streamlines (left panels), and magnitude of the vorticity ω (right panels) for the real gravel topography configuration.
Water 12 00690 g010
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 .
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 .
Water 12 00690 g011
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 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.
Water 12 00690 g012
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.
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.
Water 12 00690 g013

Share and Cite

MDPI and ACS Style

Tavelli, M.; Piccolroaz, S.; Stradiotti, G.; Pisaturo, G.R.; Righetti, M. 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. Water 2020, 12, 690. https://doi.org/10.3390/w12030690

AMA Style

Tavelli M, Piccolroaz S, Stradiotti G, Pisaturo GR, Righetti M. 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. Water. 2020; 12(3):690. https://doi.org/10.3390/w12030690

Chicago/Turabian Style

Tavelli, Maurizio, Sebastiano Piccolroaz, Giulia Stradiotti, Giuseppe Roberto Pisaturo, and Maurizio Righetti. 2020. "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" Water 12, no. 3: 690. https://doi.org/10.3390/w12030690

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop