Extending the Shields criterion to erosion of weakly cemented granular soils

This contribution tackles the issue of incipient conditions for initiation of erosion by a fluid flow at the surface of cohesive materials. To this end, a typical assessment procedure consists of subjecting a soil sample to progressive hydrodynamic stresses induced by a submerged impinging jet flow whose injection velocity is gradually increased. This paper presents the results of an extensive use of this protocol both in experiments and numerical simulations, the latter being based on a coupled DEM and LBM approach. Here we consider the specific case of weakly cemented soils, either made experimentally of glass beads bonded by solid bridges or modelled numerically by a solid bond rheology with a parabolic yield condition involving the micromechanical traction, shearing and bending of the bonds. The results show that, as expected, the hydrodynamic stress for erosion onset substantially increases with solid cohesion as compared to cohesionless cases but can, however, be satisfactorily predicted by a simple extension of the usual Shields criterion that only applies for cohesion-less granular sediments. This extension includes a cohesion number, the granular Bond number, with a simple definition based on tensile yield values. 1 Context and motivation The ability to better understand and correctly predict sediment erosion and transport is of paramount interest owing to the large number of related practical situations in nature or industry. Soil resistance to erosion by a surface fluid flow is indeed an old research topic whose foundation stone was laid almost one century ago by the pioneering work of A. F. Shields [1], mainly focused on the onset of erosion in granular sediments. However, despite a continuously growing number of contributions on the subject, there is still a lack of knowledge on the fundamental issue of the critical flow conditions for erosion of natural soils. Many of these difficulties have undoubtedly to do with the complexity of the phenomenon, both in terms of the hydrodynamic conditions (turbulent or transitional regimes, stress fluctuations in time and space) and of the sediment's nature (particle size distribution, particle shapes, internal stresses due to friction, adhesion and/or electrostatic forces), especially since real soils often present a certain degree of macroscopic cohesion. For strictly non-cohesive materials, erosion can be simply considered as a grain-by-grain process that initiates as soon as the fluid flow stress exceeds both the particle weight and the frictional forces at the sediment’s surface. This condition can be quantified by a critical value Sh0 ∗ of the so-called Shields number Sh0, which is the dimensionless ratio between the shear stress exerted by the fluid flow over the sediment and the buoyant weight of a single particle. Values of Sh0 ∗ range * Corresponding author: pierre.philippe@inrae.fr over almost one order of magnitude depending on the particle shear Reynolds number Reτ and delineate the so-called Shields curve. Several empirical explicit formulations exist for this curve, as for instance that given by Guo [2] that we will use here. In the case of cohesive materials (benthic sediments, clayey soils, cemented calcareous sands, etc), the complexity of the problem increases due to the additional contribution of internal attractive forces between particles. Substantial collective processes during erosion are thus involved and essentially invalidate the previous Shields approach. Nonetheless, some efforts to extend it can be found in the literature [3-7], mostly based on the definition of a second dimensionless number that compares a typical internal cohesive force to the sediment buoyant weight and is often denoted Bog for granular Bond number [8-9]. The objectives of the present study are multiple: (i) to extent the applicability of the Shields criterion to weakly cemented sediments; (ii) to compare real experimental results with their numerical counterparts from a 2D DEM-LBM modelling; (iii) to evaluate a characteristic cohesive stress and thus define a granular Bond number at both contact and sample scales. To this end, the experiments and the numerical approach with a common protocol are both described in section 2 while section 3 presents the results, including a comparison between experiments and simulations. To conclude we propose an extended and more general formulation for the critical Shields number involving the granular Bond number but still restricted to the present specific © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). EPJ Web of Conferences 249, 08009 (2021) https://doi.org/10.1051/epjconf/202124908009 Powders and Grains 2021 A video is available at https://doi.org/10.48448/qgg0-3y46 cohesion generated by weak cemented bonds within a


Context and motivation
The ability to better understand and correctly predict sediment erosion and transport is of paramount interest owing to the large number of related practical situations in nature or industry. Soil resistance to erosion by a surface fluid flow is indeed an old research topic whose foundation stone was laid almost one century ago by the pioneering work of A. F. Shields [1], mainly focused on the onset of erosion in granular sediments. However, despite a continuously growing number of contributions on the subject, there is still a lack of knowledge on the fundamental issue of the critical flow conditions for erosion of natural soils. Many of these difficulties have undoubtedly to do with the complexity of the phenomenon, both in terms of the hydrodynamic conditions (turbulent or transitional regimes, stress fluctuations in time and space) and of the sediment's nature (particle size distribution, particle shapes, internal stresses due to friction, adhesion and/or electrostatic forces), especially since real soils often present a certain degree of macroscopic cohesion.
For strictly non-cohesive materials, erosion can be simply considered as a grain-by-grain process that initiates as soon as the fluid flow stress exceeds both the particle weight and the frictional forces at the sediment's surface. This condition can be quantified by a critical value ℎ 0 * of the so-called Shields number ℎ 0 , which is the dimensionless ratio between the shear stress exerted by the fluid flow over the sediment and the buoyant weight of a single particle. Values of ℎ 0 * range * Corresponding author: pierre.philippe@inrae.fr over almost one order of magnitude depending on the particle shear Reynolds number and delineate the so-called Shields curve. Several empirical explicit formulations exist for this curve, as for instance that given by Guo [2] that we will use here.
In the case of cohesive materials (benthic sediments, clayey soils, cemented calcareous sands, etc), the complexity of the problem increases due to the additional contribution of internal attractive forces between particles. Substantial collective processes during erosion are thus involved and essentially invalidate the previous Shields approach. Nonetheless, some efforts to extend it can be found in the literature [3][4][5][6][7], mostly based on the definition of a second dimensionless number that compares a typical internal cohesive force to the sediment buoyant weight and is often denoted for granular Bond number [8][9]. The objectives of the present study are multiple: (i) to extent the applicability of the Shields criterion to weakly cemented sediments; (ii) to compare real experimental results with their numerical counterparts from a 2D DEM-LBM modelling; (iii) to evaluate a characteristic cohesive stress and thus define a granular Bond number at both contact and sample scales. To this end, the experiments and the numerical approach with a common protocol are both described in section 2 while section 3 presents the results, including a comparison between experiments and simulations. To conclude we propose an extended and more general formulation for the critical Shields number involving the granular Bond number but still restricted to the present specific cohesion generated by weak cemented bonds within a granular material.

Methodology
To increase the robustness of the correlations developed here we have conducted both experiments and numerical simulations with a similar fluid flow configuration and erosion test procedure as detailed below. An important issue thereby concerns the proper quantification of a characteristic value for the cohesive strength of each tested sample, either at contact-scale or at sample-scale, both quantities being theoretically related by the classical homogenisation law proposed by Rumpf [10].

Experiments with artificial materials
The samples used for the experiments are artificially bonded granular materials made of spherical glass beads with solid bridges out of polyurethane resin. As detailed in [11], the resin is initially liquid and diluted at varying concentrations in water before being mixed with the glass beads. The final samples to be tested are obtained after complete drying of the resin and transformation of the liquid bridges into solid bonds (duration of one day).
A first set of experiments, denoted E1, was carried out with beads of diameter 3.00 ± 0.02 mm made of borosilicate glass (density = 2230 kg. m −3 ). This particular choice of transparent beads makes it possible to implement the so-called refractive index matching technique (RIM) based on the immersion of the beads in a liquid that shares the same refractive index, in this case ≈ 1.472. The immersion liquid used here is an oil mixture whose density is = 846 ± 5 kg. m −3 and whose dynamic viscosity is = 28 ± 2 cP. The RIM technique allows for a better visualisation, especially during scour development once erosion is initiated [11], but this is out of the scope of the present contribution.
Two further experimental series were implemented for samples made of silicate glass beads (density = 2495 kg. m −3 ) immersed in the same liquid but with no RIM. Two bead diameters were used in the ranges 2.85 − 3.30 mm and 0.75 − 1.00 mm, denoted E2 and E3 respectively.
A key issue about these artificial materials concerns the experimental determination of their cohesive strength. To this end, direct tensile tests were carried out at the single-contact scale for all the 3mm-bead samples. After some averaging, a mean microscopic yield force can be evaluated. As such tests were not feasible for the too small 1mm-beads, a series of macroscopic tensile tests was instead performed at the sample scale to quantify in this case a macroscopic yield stress . Additional details about the different set-ups are provided in [11].

Numerical simulations with DEM-LBM
Our numerical investigation at micro-scale relies on a combination of the Discrete Elements Method (DEM) and the Lattice Boltzmann Method (LBM). The DEM is often employed for the mechanical description of assemblies of solid particles mainly interacting by contacts through friction and collisions, while the LBM can accurately simulate complex fluid flows, including interstitial flows within porous media, based on an implicit resolution of the discrete form of the Boltzmann equation that can retrieve the incompressible Navier-Stokes equation at small Mach numbers. The 2D simulation considered involves circular particles with a mean diameter and a uniform size distribution between 0.8 and 1.2 . Here the mean diameters are 1 = 2 mm, 2 = 3 mm, and 3 = 5 mm. Three values of viscosities are also used: 1 = 40 0 , 2 = 50 0 , and 3 = 100 0 , with 0 the kinematic viscosity of water. We have selected 6 different sets of parameters: S1 ( 1 , 2 ), S2 ( 2 , 2 ), S3 ( 2 , 1 ), S4 ( 2 , 3 ), S5 ( 3 , 1 ), and S6 ( 1 , 1 ). Further details can be found in [12][13].
In this model, we use a specific solid bond contact model at each particle-particle contact as proposed by Delenne and co-workers [14]. It is an elastoplastic bond model associated with three local degrees of freedom at contact points (normal, tangential and rotational) and involving yield values for the pure normal force ( ), shear force ( ), and bending moment ( ). Supported by some experimental results [14], a paraboloidal yield surface in the space of interactions is here assumed to account for bond failure by mixed solicitations. Once the bond is broken, the contact becomes purely frictional. Finally, adopting the ratios between the single-mode thresholds measured in [14], we can relate all yield parameters to a unique cohesion force defined by: (1)

Protocol for erosion onset measurement
We have used similar experimental and numerical protocols, that consider a submerged jet interacting with the surface of the sample being tested. This jet erosion test is commonly used to quantify the resistance to erosion of natural soils, either in the lab or in the field. Here, the mean injection velocity of the jet is gradually increased and the erosion onset is identified with the first grain motion. Then, direct observations during the test, together with the time evolution of the number of eroded grains in the simulations, enable an accurate measurement of the critical injection velocity * , especially for the numerical data whose error bars are smaller the symbol size in the forthcoming graphs.
Since the Shields criterion requires the determination of the critical shear-stress * exerted at the sample's surface, a relation is needed to deduce * from * . Most previous studies on turbulent impinging jet erosion have relied on empirical expressions derived from experiments of solid wall impingements or even have neglected the presence of the soil surface to permit directly the use of the well-known self-similar free jet and its analytical expression, either in 2D or 3D. In a recent study based on the same 2D LBM-DEM modelling in laminar flow regime [15], we have highlighted the adequacy of the free jet models to quantify the maximal shear-stress at the soil surface pointing out that: (i) the maximal horizontal velocity at the vicinity of the soil surface is directly proportional to the theoretical vertical velocity of the free jet model at the same distance from jet's injection; (ii) the maximal shear-stress at the soil surface is given by a simple inertial expression owing that a Blasius-like friction coefficient is added, i.e. a friction coefficient inversely proportional to the square root of the jet Reynolds number as observed for a laminar flow along a flat wall. The same type of friction term is consequently assumed for our 3D experiments and the empirical proportionality coefficient can be fixed by comparison with previous cohesionless results [11,16]. Then, the bed shear-stress at erosion onset * for both our numerical simulations and our experiments is calculated from the critical injection velocity * via these 2D and 3D expressions, respectively.

Inadequacy of the Shields curve
Data from several experimental tests and numerical simulations of erosion tests are represented in the two Shields diagrams of Figure 1 and Figure 2. In these graphs, the critical Shields number ℎ 0 * = * (∆ ) ⁄ is plotted versus the critical particle shear Reynolds number * = √ * ⁄ ⁄ , where ∆ stands for the difference between the grain and the fluid densities, and is the mean diameter of the particles. symbols) and the extended * (red closed symbols) Shields numbers as a function of the critical particle shear Reynolds number * deduced from the experiments [11]. The solid line stands for the so-called Shields curve [2].
The expression of the Shields number used here is the original one, developed for granular sediments, and, as can be obviously noticed on the two graphs, a systematic and increasing deviation from the Shields curve (here given by the implicit expression of Guo [2]) is found when * increases. The standard onset description through the Shields curve therefore appears misleading for our weakly cemented materials, especially when the adhesion of the solid bonds increases. This finding thus calls for a new criterion that incorporates adhesion at contacts, in addition to weight and friction, to possibly account for our cemented granular materials.

Extension of Shields criterion approach
A straightforward way to extend the Shields number definition is simply to add the contribution from the solid bond adhesion to the other stresses that are considered to resist erosion. These include the buoyant weight and friction, the latter being just proportional to the former. This extension requires a characteristic value in the form of a cohesive stress . When a yield tensile stress is measured at the sample scale, one can directly identify with . On the contrary, the passage from microscopic forces, either obtained experimentally for large enough beads or defined as the characteristic yield load within the numerical contact bond model, to requires an homogenisation relation. Relying on the expression by Rumpf [10] and introducing realistic values for the coordination number in a bead pack [11], the following relation can be proposed: where the microscopic force stands both for (experiments with 3mm-beads) and (numerical simulations).
Then, the extended expression of the Shields number reads: is the granular Bond number which compares cohesion to buoyant weight.
In this expression we have also introduced an empirical coefficient , expected to be of order 1, that arises from pre-factors and is mainly ruled by geometry. Based on the values obtained for the two critical Shields numbers ℎ 0 * and ℎ * , the coefficient can be determined quantitatively by linear regression from: ℎ * ℎ 0 * = 1 + where the implicit expression of Guo [2] is used for the Shields curve given by ℎ 0 * as a function of * . As shown in Figure 3, linear relations are indeed in correct agreement with our data and provide two different but unexpectedly rather close values for : 3 = 2.26 ± 0.27 and 2 = 1.82 ± 0.03 for the experimental and numerical results, respectively. Using these values for , the critical values of the extended Shields number ℎ * can now be calculated and match well with the Shields curve, as shown in Figures 1 and  2. This finding demonstrates the benefits of using ℎ * , which improves significantly the predictions of the usual expression ℎ 0 * of the Shields number while reconciling the data for cohesive and cohesion-less materials.

Conclusion
This study presents a series of physical jet erosion tests and complementary data from 2D micromechanical simulations and shows that, unsurprisingly, the traditional formulation of the Shields number used for cohesion-less sediments is inadequate to account for the onset of erosion in weakly cemented soils. We therefore propose an extension of the usual Shields number that, in addition to the buoyant weight and friction, also includes adhesion by solid bonds at contacts through either a yield tensile force at the micro-scale or a yield tensile stress at the sample-scale. Based on Rumpf's theory [10], a unique definition can then be used for the granular Bond number , comparing cohesion to buoyant weight, leading this way to the expression given in Equation (3) for the extended Shields number ℎ. The comparison with the present experimental and numerical data demonstrates a rather satisfactory agreement with the usual Shields curve, with respective values of coefficient approximately equal to 2.3 in 3D (experiments) and 1.8 in 2D (numerical simulations). Note however that other factors than dimensionality could produce differences between experiments and simulations. This work thereby offers an extended framework that is valid for both cohesionless (for ≪ 1) and cemented materials (for > 1), at least for the weakly cemented samples studied here, with < 300.