# Extending the Shields criterion to erosion of weakly cemented granular soils

Pierre Philippe, Zeyd Benseghier, Florian Brunier-Coulin, Li-Hua Luu, Pablo Cuéllar, Stéphane Bonelli

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Philippe, Zeyd Benseghier, Florian Brunier-Coulin, Li-Hua Luu, Pablo Cuéllar, et al. Extending the Shields criterion to erosion of weakly cemented granular soils. EPJ Web of Conferences, 2021, 249, pp.08009. 10.1051/epjconf/202124908009. hal-04116508

HAL Id: hal-04116508

https://hal.inrae.fr/hal-04116508

Submitted on 4 Jun 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extending the Shields criterion to erosion of weakly cemented granular soils

Pierre Philippe<fnref n="1" />*, Zeyd Benseghier<fnref n="1" />, Florian Brunier-Coulin<fnref n="1" />,<fnref n="2" />, Li-Hua Luu<fnref n="1" />, Pablo Cuéllar<fnref n="3" />, and Stéphane Bonelli<fnref n="1" />

**Abstract.** The present 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 in 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 cohesion-less 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 erosion onset for granular sediments. However, despite a continuously growing number of contributions on the subject, there is still a lack of knowledge on such a fundamental issue as 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 $Sh_0^*$ of the so-called Shields number $Sh_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 $Sh_0^*$ range

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 the one given by Guo [2] that we will use in the following.

By considering now the 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 $Bo_g$ 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

<footnote n="1">¹RECOVER, UMR 1467, INRAE and Aix Marseille Univ, Aix-en-Provence, France</footnote>

<footnote n="2">²Irphé, UMR 7342, CNRS and Aix-Marseille Univ, Marseille, France</footnote>

<footnote n="3">³BAM, Federal Institute for Materials Research and Testing, Berlin, Germany</footnote>

<footnote>* Corresponding author: pierre.philippe@inrae.fr</footnote>
cohesion generated by weak cemented bonds within a granular material.

# 2 Methodology

As explained just before, the main goal of this contribution is to test and improve the semi-empirical approach based on the Shields criterion when surface erosion is induced by a controlled fluid flow for artificial cohesive soils, in this case weakly cemented materials. To increase the robustness of the methodology 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].

## 2.1 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.

A first set of experiments was carried out with beads of diameter $3.00 \pm 0.02$ mm made of borosilicate glass (density $\rho_b = 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 $n_r \approx 1.472$. The immersion liquid used here is an oil mixture whose density is $\rho_l = 846 \pm 5$ kg.m$^{-3}$ and whose dynamic viscosity is $\mu_l = 28 \pm 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 $\rho_s = 2495$ kg.m$^{-3}$) immersed in the same liquid but with no RIM. Two bead diameters were used in the ranges $0.75 - 1.00$ mm and $2.85 - 3.30$ mm, 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 $F_m^Y$ 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 $\tau_M^Y$. Additional details about the different set-ups are provided in [11].

## 2.2 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 present 2D simulation involves circular particles with a mean diameter $d$ and a uniform size distribution between $0.8d$ and $1.2d$. Further details can be found in [12-13].

In this model, we use a specific solid bond rheology at each particle-particle contact as proposed by Delenne and co-workers [14]. It is an elastoplastic bond rheology associated with three local degrees of freedom at contact points (normal, tangential and rotational) and involving yield values for the pure normal force ($F_n^Y$), shear force ($F_t^Y$), and bending moment ($M_b^Y$). 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 $F_c$ defined by:

$$F_c = F_n^Y = 2F_t^Y = 4 \frac{M_b^Y}{d} \quad (1)$$

## 2.3 Protocol for erosion onset measurement

We have used similar experimental and numerical protocols, consisting in the hydrodynamic solicitation at the surface of a tested sample by a submerged impinging jet. Such a 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 $u_j^*$.

Since the Shields criterion requires the determination of the critical shear-stress $\tau_f^*$ exerted at the sample's surface, a relation is needed to deduce $\tau_f^*$ from $u_j^*$. 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 cohesion-less results [11, 16]. Then, the bed shear-stress at erosion onset $\tau_f^*$ for both our numerical simulations and our experiments is calculated from the critical injection velocity $u_j^*$ via these 2D and 3D expressions, respectively.

# 3 Results and discussion

## 3.1 Inadequacy of the Shields curve

The data of several experimental and numerical series of such erosion tests are represented in the two Shields diagrams of Figure 1 and Figure 2. In these graphs, the critical Shields number $Sh_0^* = \tau_f^*/(\Delta\rho gd)$ is plotted versus the critical particle shear Reynolds number $Re_\tau^* = \rho_l d\sqrt{\tau_f^*/\rho_l}/\mu_l$, where $\Delta\rho$ stands for the difference of density between the grains and the fluid, and $d$ is the mean diameter of the particles.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Onset values of both the usual $Sh_0^*$ (black open symbols) and the extended $Sh^*$ (red closed symbols) Shields numbers as a function of the critical particle shear Reynolds number $Re_\tau^*$ deduced from the experiments. The solid line stands for the so-called Shields curve [2].</figcaption></figure>

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 diagrams, a systematic and increasing deviation from the Shields curve (here given by the implicit expression of Guo [2]) is found when $Re_\tau^*$ 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.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Onset values of both the usual $Sh_0^*$ (black open symbols) and the extended $Sh^*$ (red closed symbols) Shields numbers as a function of the critical particle shear Reynolds number $Re_\tau^*$ obtained with the 2D numerical simulations. The solid line stands for the so-called Shields curve [2].</figcaption></figure>

## 3.2 Extension of Shields criterion approach

A straightforward way to extend the Shields number definition is simply to add a contribution for solid bond adhesion in the inventory of all resistant stresses against erosion, which already includes 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 $\sigma_c$. When a yield tensile stress $\tau_M^Y$ is measured at the sample scale, one can directly identify $\sigma_c$ with $\tau_M^Y$. On the contrary, the passage from microscopic forces, either $F_m^Y$ obtained experimentally for large enough beads or $F_c$ defined as the characteristic yield load within the numerical contact bond model, to $\sigma_c$ 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:

$$ \sigma_c = 1.25 F/d^2. \qquad (2) $$

where the microscopic force $F$ stands both for $F_m^Y$ (experiments with 3mm-beads) and $F_c$ (numerical simulations).

This way, the extended expression of the Shields number reads:

$$ Sh = \frac{\tau_f}{\Delta\rho gd + \alpha \sigma_c} = \frac{Sh_0}{1 + \alpha Bo_g} \qquad (3) $$

where $Bo_g = \sigma_c/(\Delta\rho gd)$ is the granular Bond number which compares cohesion to buoyant weight.

In this expression we have also introduced an empirical coefficient $\alpha$, expectedly 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 $Sh_0^*$ and $Sh^*$, the coefficient $\alpha$ can be determined quantitatively by linear regression from:

$$ \frac{Sh^*}{Sh_0^*} = 1 + \alpha Bo_g \qquad (4) $$
where the implicit expression of Guo [2] is used for the
Shields curve given by $Sh_0^*$ as a function of $Re_τ^*$.

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 $\alpha$: $\alpha_{3D} = 2.26 \pm 0.27$ and $\alpha_{2D} = 1.82 \pm 0.03$ for the experimental and numerical results, respectively. Using these values for $\alpha$, the critical values of the extended Shields number $Sh^*$ 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 $Sh^*$, which improves significantly the predictions of the usual expression $Sh_0^*$ of the Shields number while reconciling the data for cohesive and cohesion-less materials.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Plot and linear regression (with the corresponding values of $\alpha$) of $(\frac{Sh^*}{Sh_0^*} - 1)$ versus the granular Bond number $Bo_g$ for (a) the 3D experiments and (b) the 2D numerical simulations.</figcaption></figure>

## 4 Conclusion

This study presents several series of jet erosion tests providing both experimental data with artificial materials and numerical data from 2D micromechanical models and shows that, unsurprisingly, the traditional formulation of the Shields number used for cohesion-less sediments is inadequate to account for the erosion onset of 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 $Bo_g$, comparing cohesion to buoyant weight, leading this way to the expression given in Equation (3) for the extended Shields number $Sh$. The comparison with the present experimental and numerical data demonstrates a rather satisfactory agreement with the usual Shields curve, with respective values of coefficient $\alpha$ approximately equal to 2.3 in 3D (experiments) and 1.8 in 2D (numerical simulations). This work thereby offers an extended framework that is valid for both cohesion-less (for $Bo_g \ll 1$) and cemented materials (for $Bo_g > 1$), at least for the weakly cemented samples studied here, with $Bo_g < 300$.

Two of the authors, FBC and ZB, are grateful for the financial support provided by the "Région Sud, Provence-Alpes-Côte d'Azur". The authors would also like to thank P. Aussillous, J.-Y. Delenne, J. Duriez, F. Lominé, and Y. Khidas for fruitful discussions. The contribution of J.-L. Sarrat as an internship student is also acknowledged

## References

1. J. M. Buffington, J. Hydraul. Eng. **125**, 4 (1999)

2. J. Guo, J. Hydraul. Eng. **123**, 7 (1997)

3. W. Lick, L. Jin, J. Gailani, J. Hydraul. Eng. **130**, 8 (2004)

4. P. Claudin, B. Andreotti B, Earth Planet. Sc. Lett. **252**, 1-2 (2006).

5. M. Righetti, C. Lucarelli, J. Geophys. Res. **112** C05039 (2007)

6. F. Ternat, P. Boyer, F. Anselmet, M. Amielh, Water Resour. Res. **44**, W11434 (2008).

7. D. Chen, Y. Wang, B. Melville, H. Huang, W. Zhang, J. Hydraul. Eng. **144**, 4 (2018)

8. S. T. Nase, W. L. Vargas, A. A. Abatan, J. J. McCarthy, Powder Technol. **116** (2001)

9. A. Anand, J. S. Curtis, C. R. Wassgren, B. C. Hancock, W. R. Ketterhagen, Chem. Eng. Sci. **64** (2009)

10. P. Pierrat, H. Caram, Powder Technol. **91** (1997)

11. F. Brunier-Coulin, P. Cuéllar, P. Philippe, Phys. Rev. Fluids **5**, 3 (2020)

12. J. Ngoma, P. Philippe, S. Bonelli, F. Radjaï, J.-Y. Delenne, Phys. Rev. E **97**, 5 (2018)

13. Z. Benseghier, P. Cuéllar, L.-H. Luu, S. Bonelli, P. Philippe, Comput. Geotechn. **120** (2020)

14. J.-Y. Delenne, M. S. El Youssoufi, F. Cherblanc, J.-C. Bénet, Int. J. Numer. Anal. Meth. Geomech. **28**, 15 (2004)

15. Z. Benseghier, P. Cuéllar, L.-H. Luu, J.-Y. Delenne, S. Bonelli, P. Philippe, J. Hydraul. Eng. **141**, 4 (2019)

16. F. Brunier-Coulin, P. Cuéllar, P. Philippe, Phys. Rev. Fluids **2**, 3 (2017)