Erosion of cohesive soil layers above underground conduits

Using a recently developed 2D numerical modelling that combines Discrete Element (DEM) and Lattice Boltzmann methods (LBM), we simulate the destabilisation by an hydraulic gradient of a cohesive granular soil clogging the top of an underground conduit. We aim to perform a multi-scale study that relates the grain scale behavior to the macroscopic erosion process. In particular, we study the influence of the flow conditions and the inter-particle contact forces intensity on the erosion kinetic.


Introduction
Location of buried cavities, mostly in urban zones, represents a highly vulnerable area with regard to geotechnical aspects. Their collapse, that leads to sinkholes, is indeed a major cause of large soil subsidences. Our research focuses more specifically on the destabilization, after floods, of weak layers above existing conduits due to underground water flows, as depicted in Figure 1a. Up to now, the understanding of the resistance against erosion of cohesive soils remains poor, being mainly limited to empirical laws based on experimental tests performed in laboratory or insitu [1,2]. In this study, we propose to implement a numerical model, increasingly used in geomechanics [3,4], that combines the Lattice Boltzman Method (LBM) and the Discrete Element Method (DEM) to efficiently describe the two-phase flow. As demonstrated in recent works, the proposed numerical model has the advantage to provide access to the local erosion behavior, at the grain scale [5][6][7]].

Numerical model
On the one hand, the fluid phase is simulated by LBM, which is an explicit 2D finite difference scheme over a discrete lattice mesh of the continuous Boltzmann equation that involves the collision and propagation of fluid particles. We employ the classical 2D lattice with 9 particle velocity vectors, the so-called D2Q9 scheme. The collision process is described by the multiple relaxation time model (MRT) to overcome the shortcomings of the basic Bhatnaggar-Gross-Krook scheme (BGK) [5,6].
On the other hand, the solid phase is represented by a cohesive frictional granular material using a 2D Discrete Element method [8,9]. To this end, we use a simple elastoplastic model characterized by the normal and shear bond   Table 1. stiffnesses, k n,b and k s,b , that includes a parabolic yield surface in the space of contact forces, as sketched in Figure 2. When a contact force reaches the yield surface, the bond is broken and the contact becomes purely frictional, via a Kelvin-Voigt relationship to model the normal force and a viscous-regularized Coulomb law for the shear force (see stiffness, damping and friction coefficients in Table 1). In the present numerical model, the normal and shear yield thresholds C n and C s are assumed to depend only on a single parameter C=C n = 2·C s , which represents the inter- Figure 2. Yield surface of cohesive bonds in the space of interaction normal and shear forces (F n , F s ), from [6] particle bond strength. The soil cohesion can then be characterized by a dimensionless cohesion number defined as the ratio of bond strength to the particle's own buoyant weight, Coh = C/(ρ s − ρ f )gV, where ρ s − ρ f is the submerged apparent density, g is the gravitational acceleration and V is the particle's volume.
To deal with the solid-fluid coupling, we use a relationship based on a momentum-exchange method [10]. Finally, the LBM calculation is implemented with a hydraulic radius smaller than the real particule radius to recover a non-zero permeability in our 2D description.
Using this combined DEM-LBM method, we thus intend to model a specific soil erosion situation. The simulation setup consists in a conduit initially clogged by an immersed granular assembly as sketched in Figure 1b. In order to impose a preferential path for the flow within the granular sample, two soil cohesions are settled. In the external layers, the bond strength is fixed at C=1000 N, whereas in the central layer we set a much weaker cohesion magnitudes from 5 to 20 N. The erosion of the soil layer is triggered by a pressure gradient. As listed in Table  1, we simulate two samples of D=6 mm diameter grains, with a same height but two different lengths L, corresponding to particle numbers N p equal to 11816 and to 3042.

Erosion regimes
The behavior of the granular layer submitted to a gradient pressure is parametrically studied. Table 1 reports the parameters implemented. By varying the conduit size ∆ defined by ∆ = l/D, the granular cohesion given by the cohesion number Coh and the inlet fluid pressure P, we observe a "blocking" stable case where no erosion occurs and an intermittent evacuation of the granular sample through the predefined internal conduit. We here propose to separate this phenomenon into three different erosion regimes, as presented in Figure 3. At the lowest magnitude of erosion, a stable cavity can marginally be formed. If increasing eroding parameters, a backward erosion mechanism arises, where the soil destabilisation front direction is opposite to the flow direction. Once the erosion front reaches the top of the sample, the detached matter is evacuated by the flow. Finally, we call plug extrusion, the regime where    the whole granular column falls all at once. It is important to note that one of the LBM limitation is that the fluid velocity has to be much smaller than the lattice speed c s , so as to remain in the incompressible Navier-Stockes behavior. In this study, calculations stop when the fluid reaches 1 m/s. Consequently, we often cannot simulate the complete evacuation from the conduit, but as examinated later, its kinetic seems to reach a constant velocity that we are able to measure. In Figure 4, we plot the three phase diagrams related to ∆, Coh and P for the two samples. We notably indicate the erosion threshold by the lines. Globally, erosion process can be triggered by a high enough conduit size, a high enough inlet pressure and a small enough cohesion. The first diagram shows that the higher Coh is, the higher ∆ is required to erode soil. As regards the eroding effect of increasing P, the two other diagrams display very sharp threshold. This observation, that we are currently studying more thoroughly, could be related to Janssen effects.
Our simulations also point out the influence of the granular sample length L on the erosion threshold. This latter is slighthy "delayed" for the smaller sample, in the sense that a higher critical ∆, a smaller critical Coh and a higher critical P are involved in the erosion process. To capture the physics underlying this result, we looked at the initial state of inter-particle cohesion, just before the first breaking of bonds. Figure 5 illustrates the normal forces between particles at this stage for both samples. We can distinguish a lower magnitude in the smallest sample specifically within the relevant central layer, which means a greater distance to the yield threshold (presented in Fig.  2). We can notice the "conic" shape of the normal force distribution, that is even more pronounced with a higher sample length.

Erosion kinetic
From simulation pictures, we propose to plot spatiotemporal diagrams to analyse the erosion kinetic. Figure  6a displays a typical analysis following 3 steps (from left to right): (i) selection of the central part of the granular assembly, where particle motion occurs ; (ii) proceeding on spatio-temporal diagrams within this area along each vertical line with a spatial resolution of 1 mm ; (iii) calculation of the standard deviation over the set of diagrams. In the  first case of backward extrusion, the time evolution of the grain's collective motion clearly presents two different kinetics (Fig. 6a). It first shows that the intiation of motion starts at the bottom of the granular sample while the top part stay static. This so-called backward mechanism is depicted by the green lines, horizontal at the top and tilted at the bottom. Second, particles move almost together with a constant velocity, represented by the two red lines having the same slope. Now let us slightly increase the conduit size to get the transition towards the purely plug extrusion. Figure 6b shows a typical spatio-temporal diagram of this regime. The red line indicates that we could consider, in a first approximation, a roughly constant evacuation of the whole block with the same velocity than during the final plug-like motion of the backward extrusion. This result therefore demonstrates a continuous transition between these two erosion processes.

Conclusion
This paper presents first results of a parametric study on numerical simulations of the destabilisation of a soil laying above an underground conduit due to a pressure gradient. We mainly discussed the influence of the inlet fluid pressure, the conduit size and the inter-particle cohesion  on the different erosion regimes. One of the important findings is the role of the granular sample length on the erosion threshold. Complementary ongoing studies will consider effects of the particle size and the granular sample height, in the aim of capturing the relevant parameters for predicting the different erosion situations. Exploiting the grain scale behavior, we also analysed both the internal cohesion state of the soil and the erosion kinetic. More generally, this numerical study confirmed the relevance of employing a DEM-LBM coupling to get insight into micromechanical features of soil erosion realted to complex natural situations such as blocked karst conduits.