Engineering Geology 260 (2019) 105249

Contents lists available at ScienceDirect

Engineering Geology
journal homepage: www.elsevier.com/locate/enggeo

Hydro-mechanical modeling of sinkhole occurrence processes in covered T

karst terrains during a flood
⁎ ⁎
Li-Hua Luua, , Gildas Nouryb, , Zeyd Benseghiera, Pierre Philippea
a
Irstea, Aix Marseille Univ, RECOVER (Risks, Ecosystems, Vulnerability, Environment and Resilience), 3275 route de Cézanne, 13182 Aix-en-Provence, France
b
BRGM (French Geological Survey), DRP (Risks and Prevention Division), 3 avenue Claude Guillemin, BP 36009, Orléans 45060, France

A R T I C LE I N FO A B S T R A C T

Keywords: Groundwater is an aggravating factor in karstic sinkhole activity as it exacerbates infiltration, percolation, soil
Cover-collapse sinkhole saturation and drainage. The exceptionally high number of ground collapses triggered during the major flood of
Flood spring 2016 in the Orléans region (France) clearly supports this assertion. In this article, we examine the role of
Karst terrain flooding in sinkhole occurrence in cohesive soil layers covering karstified limestone rock. An innovative hydro-
Cohesive soil erosion
mechanical model is applied to simulated field scenarios. Our numerical simulations combine the Discrete
Hydro-mechanical modeling
Element Method (DEM) to model the solid phase with the Lattice Boltzmann Method (LBM) for the fluid phase.
This coupled numerical method allows us to explore the micromechanical features of internal soil erosion in a
flood situation. Three processes, consistent with field observations, are simulated and studied through phase
diagrams: the formation of a stable cavity within the cover material, the upward propagation of a cavity leading
to a dropout sinkhole, and the downward discharge of the granular media, called the subsidence sinkhole. In
particular, we perform a parametric analysis of the dropout sinkhole that gives an estimate of the collapse width.
In the first approximation, the characteristic length is shown to increase linearly with cover thickness, regardless
of the other main parameters (soil cohesion, hydraulic head, system geometry). Finally, we present an ex-
ploratory experimental study, using sand with artificial cohesion that successfully reproduces different erosion
regimes predicted by our numerical simulations.

1. Introduction limestone cavities via water circulation (Beck, 2012). Moreover, karst
cover soil is particularly hazardous in the case of intense rainfall and
Karstic landscapes cover around 15–20 % of the Earth's surface and floods (Brinkmann et al., 2008; Gutiérrez et al., 2014; Hyatt and Jacobs,
are a matter of concern for land-use, especially for sinkhole hazards 1996). In France, the meteorological event that occurred in spring 2016
(Ford and Williams, 2013). Cover by limestone karst is one of the most is a new example of the impact of intense rainfalls and floods on
critical environments (Sowers, 1996; Waltham, 2008), mostly located sinkhole occurrence (Noury et al., 2018). The lessons learned from this
in easy-to-build flat areas (such as alluvial plain). The subterranean event enabled the French Geological Survey (BRGM) and the National
configuration is often inaccessible (filled or underwater entries) while Environment and Agriculture Research Institute (IRSTEA) to improve
the non-karstic lithology of the cover seems “stable”. Large urbanized insight into internal erosion processes through innovative numerical
areas have been built on this kind of karst cover soil, the most promi- modeling.
nent examples being Florida in the USA, the Gauteng Province in South This article originated from this case study, thus it focuses on
Africa, and the Orléans area in France. Hazard and risk assessment for sinkhole formation during floods. This process is still poorly under-
the safety of people and property appears difficult and remains limited stood: some experiments have succeeded in reproducing simple con-
to the locations of recent sinkholes and the use of archives to localize figurations (Perez et al., 2017) but numerical models are rare or do not
old refilled sinkholes. The Florida Office of Insurance Regulation esti- take into account the complexity of the influence of water circulation.
mated an average economic loss of 280 million per year due to sink- Few studies address the complexity of cavity development within the
holes during the period 2006 to 2010 (Kuniansky et al., 2016). cover deposit. The main approach adopts continuum mechanics
The main process leading to sinkholes in karst cover is soil piping, schemes (often using FLAC software) to simulate the geomaterial
defined as down-washing of the material cover into the network of system using constitutive models, typically incorporating visco-elasto-

⁎
Corresponding authors.
E-mail addresses: li-hua.luu@irstea.fr (L.-H. Luu), g.noury@brgm.fr (G. Noury).

https://doi.org/10.1016/j.enggeo.2019.105249
Received 31 January 2019; Received in revised form 30 July 2019; Accepted 3 August 2019
Available online 13 August 2019
0013-7952/ © 2019 Elsevier B.V. All rights reserved.
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

plastic rheology (Baryakh and Fedoseev, 2011; Shalev and Lyakhovsky, rainfalls, about 140 mm of water was recorded over 5 days, the
2012). With the same type of numerical method, recent studies have equivalent of three normal rainy months. 2) The overflow of an artifi-
benefitted from adding hydro-mechanical coupling to take into account cial canal spread the flood over a 1-km2 area. 3) The Loire floodplain
water leakage and drawdown (Rawal et al., 2017; Tao et al., 2015). was flooded by 1 to 2 m of water, covering a recently built neighbor-
Although this kind of analysis allows modeling the general failure of an hood, the water-treatment plant and several fields for 10 days. 4) This
opening void (Mohr-Coulomb envelope), it cannot predict local failures exceptional situation led to the massive infiltration of water. The small
such as block dislocation. town of Chécy was particularly hard hit. There, twelve sinkholes
All fluid-soil interactions in cover-sinkholes are related to the formed immediately or a few days after the meteorological event
widespread phenomenon of soil erosion. In particular, the expansion of (Fig. 1). Considering the geological terrains affected and whose main
underground cavities is very similar to the backward erosion piping components are presented in Fig. 2, most of these sinkholes were as-
observed in the foundations of hydraulic works (Bonelli, 2013). In this sociated with karst collapses.
field of application, studies are at a more advanced stage, notably in
experimental research (Sellmeijer et al., 2011; Wilson et al., 2015). 2.2. Detailed analysis of flood-induced sinkholes in the Chécy area
With regard to numerical investigation, the last few years have seen
substantial progress in the micromechanics of geomaterials which allow The frequency of flood-condition sinkhole occurrence (spatial and
the discrete description of soils, especially with the Discrete Element temporal) is unusually high: 12 sinkholes in 1 km2 during the 10-day
Method (DEM) developed by Cundall and Strack (1979). DEM, based on flood correspond to a sinkhole occurrence rate of 450 events per square
the realistic mechanical modeling of a large number of individual kilometer per year. The sinkhole occurrence rate evaluated in normal
particles, has become one of the most increasingly used methods for the condition (0.02 to 0.03 per square kilometer per year), is thus multi-
numerical simulation of soils in geomechanics (Radjaï and Dubois, plied by a factor from 15,000 to 22,000 because of flooding.
2011). To this end, coupling schemes have been specifically developed Fig. 3 shows pictures from BRGM field surveys after floods. We can
to obtain fluid flows at the pore scale (Cuellar et al., 2015; Lominé distinguish two different “post-mortem” shapes of sinkholes: the hour-
et al., 2013; Luu et al., 2017; Ngoma et al., 2018; Tran et al., 2017). glass shape and the inverted bowl/ teardrop shape usually associated
In this paper, we study sinkhole occurrence in covered terrains, with cover material that is non-cohesive and cohesive, respectively
using both field observations after floods and micromechanical mod- (Waltham, 2008). The soil of Val d'Orléans is mainly a slightly cohesive
eling of the phenomenon at the scale of fluid-particle interaction. clayey-silty sand and both shapes are usually observed in normal con-
Section 2 reports on the case study of the Orléans region which was ditions. In the area studied, most of the twelve sinkholes triggered by
subject to massive rainfall in 2016. An unusual number of sinkholes the 2016 flood had an hourglass shape.
surveyed in one of the flooded areas and hypothetical scenarios are Fig. 4 shows the mean diameter of the twelve sinkholes as a function
presented. Section 3 is dedicated to our numerical modeling of cohesive of floodwater level. They all have a similar size, except two singular
soil layer destabilization triggered by hydraulic load. After presenting bigger sinkholes. One of the latter is located in an intense water-flow
the coupled numerical method that combines the Lattice Boltzmann concentration area and the other is assumed to have resulted from the
Method (LBM) for describing the fluid phase and the Discrete Element coalescence of two neighboring sinkholes, explaining their larger size.
Method (DEM) for the solid phase, we present typical simulations that So, apart from these two cases, we can therefore consider that the
determine the phenomenology of the sinkhole formations. Different floodwater level does not seem to substantially influence sinkhole
analyses from a parametric study are proposed, notably based on phase diameter, which is between 0.5 m and 3 m. By contrast, the mean dia-
diagrams for the direct estimation of collapse. Then, we propose an meter of holes in normal conditions is between 1 and 2 m. Mention can
experimental method using model granular media able to qualitatively also be made of a marginal case: the biggest recent sinkhole was 16 m in
reproduce the processes observed in simulations. We finally end the diameter and 7.5 m in depth. It occurred in 2010 without flood at Saint-
paper with a broad conclusion and perspectives. Pryvé-Saint-Mesmin and destroyed a house, fortunately without ca-
sualties.
2. Case study: the 2016 meteorological crisis on sinkholes activity
in the Orléans area (France)
2.3. Field scenarios of sinkhole processes
2.1. General context
The flood-induced downward water flow within the soil is assumed
to have increased the suffosion phenomenon,1 which consists in the
Orléans is a city in the Loire Valley located 130 km south of Paris
down-washing of alluvium by water towards cavernous bedrock with
(Fig. 1). The geological bedrock of the Orléans area is constituted by the
relatively small and poorly connected voids (not necessarily well-de-
Beauce limestone (50 to 90 m thick Tertiary lacustrine limestone)
veloped and integrated in karstic cave networks). From this elementary
(Lorain, 1973). The limestone of the Loire floodplain is mantled by
mechanism of soil migration, we hypothesize two different erosion
Quaternary alluvium (5 to 15 m thick) and the groundwater level is
processes leading to the two different final sinkhole shapes presented
close to the ground surface (5 m deep). These areas are known to be
above. Fig. 5 represents a cohesive soil layer over karstic conduits in the
underlain by an extensive karstic cave network (Perrin et al., 2015).
flood-condition. The first scenario is similar to the dropout sinkhole,
However, our knowledge of network configuration is poor due to lim-
commonly reported in the classical non-submerged case, where a cavity
ited speleological investigations (few entrances, water-filled karst)
expands from the underground up to the collapse of the roof, leading to
(Moreau, 2002). The current activity in the Beauce limestone karst is
an inverted bowl sinkhole shape. Geotechnical results commonly used
attested by sinking streams, spring locations, caves found in bore holes,
by engineers for the construction of foundations in karstic terrain, show
and the occurrence of sinkholes. Since 1903, assessments have reported
that the potential collapse diameter is 2/3 the height of the overlying
approximately 640 sinkholes in the 170-km Loire floodplain around
soil layer (Sowers, 1996). In a flood situation, this prediction could
Orléans, usually called “Val d'Orléans” (Perrin et al., 2015). The latest
change. Notably, the water load could trigger the final collapse of the
report mentions a rate of 3 to 4 per year. Mean sinkhole frequency is
thin residual layer. The second scenario is similar to that of the
therefore approximately 0.02 to 0.03 occurrences per square kilometer
per year.
At the beginning of June 2016, the northern part of France, in- 1
Suffosion is distinguished from suffusion, which is an internal erosion me-
cluding the Orléans area, endured a major meteorological crisis. Fig. 2 chanism that implies the detachment of fine particles in the pore space of a
provides an illustration of the cascading effects: 1) During the massive matrix composed of larger particles.

2
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 1. Map of general context of the Orléans (Chécy area) indicating sinkholes triggered by the 2016 flood.

Fig. 2. Geological cross section through the Loire floodplain in Chécy area and sketch of the cascading effects caused by the 2016 flood (see text).

subsidence sinkhole, in which soil destabilization starts at the surface 3. Numerical modeling
and results in an hourglass shape. It is noteworthy that this process is
usually attributed to cohesionless soil (De Waele et al., 2011; Waltham, This section presents a numerical study of sinkhole formation trig-
2008). However, the following modeling work will demonstrate that gered by floods, according to previous scenarios deduced from the in-
this unexpected scenario (for cohesive soil) is quite plausible. Numer- situ observations presented in Fig. 5. To tackle this problem from a
ical simulations will show that exceptional vertical water flow might be physical perspective, we propose to bridge the gap between the mi-
responsible for intense erosion according to this second process. cromechanical phenomena at the grain scale and the macroscopic
process of the cover-collapse sinkhole, by numerically investigating the
interactions between a fluid phase and the solid phase for a large as-
sembly of bonded particles. The main objective is to explore the

3
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 3. “Post-mortem” shapes of flood-induced sinkholes in Chécy area.

15 is used here to implement the mechanical behavior of a 2D assembly of

circular particles whose trajectories are described by Newton's equa-
Singular sinkholes
tions of motion. The particles interact through contact forces, each of
intense water
coalescence of concentration area them being decomposed into a normal force Fn and a tangential force Ft
Sinkhole diameter (m)

2 sinkholes (see Fig. 6). These two components are given by a viscoelastic Kelvin-
10
Voigt model and a viscous-regularized Coulomb law, respectively. Also,
the interaction moment M is defined by the tangential force with the
particles' radii as lever arms. To account for interparticle cohesion, we
consider a parabolic yield volume in the space of contact forces and
5 moment, as developed by (Delenne et al., 2004). In the present nu-
merical model, traction, shear and bending yield thresholds are as-
sumed to depend only on a unique force C which represents the mean
interparticle bond strength. Finally, the model is enriched by an addi-
0 tional damage model adapted from Silvani et al. (2009), which pre-
0 0.5 1 1.5 2 2.5 scribes a progressive degradation of cohesive strength for subcritical
Floodwater level (m) stresses, namely those contained within the yield volume. For more
details, the reader can refer to a previous numerical study performed
Fig. 4. Influence of floodwater level on sinkhole diameter analysis in the Chécy
with this DEM cohesion model (Cuellar et al., 2015).
area.

different erosion regimes by varying the soil properties, including co- 3.1.2. Fluid phase modeling: Lattice Bolzmann Method (LBM)
hesion, and hydrodynamic conditions through a parametric study. The LBM method, based on Boltzman's kinetic theory of gases, de-
scribes the statistical behavior of a large population of fictitious fluid
3.1. Micromechanical approach particles moving randomly. The fluid flow is simulated by solving the
discretized Boltzmann equation for the probability density function
Advanced research in geomechanics demonstrated that a subtle f (→
x ,→
c , t ) to find a particle at a position → →
x with the lattice speed c at
difference in microstructure could lead to significant changes at the the given time t. The fluid dynamics occurs as the result of propagation
macroscopic scale, notably by triggering hydro-mechanical instabilities and collision of the fluid particles over a discrete regular lattice. To
such as gravitational slope failure (Iverson, 2000) and liquefaction compute this process, the present code uses a Two-Relaxation-Time
(Benahmed et al., 2004). To study soil erosion involved in sinkhole model (TRT) (Talon et al., 2012) and a classical D2Q9 scheme that
formation, we used a 2D numerical model that combines the Lattice involves a finite 2D lattice grid with 9 directions at each grid point
Boltzmann Method (LBM) for describing the fluid phase and the Dis- (including that staying at the current position, see Fig. 6) (Lallemand
crete Elements Method (DEM) for the solid phase. An increasing and Luo, 2000; Yu et al., 2003). For a given kinematic viscosity ν and
number of numerical studies rely on this coupling to treat multi-scale fluid density ρf and an adapted choice of LBM relaxation times, a
geo-mechanical problems (Cuellar et al., 2015; Lominé et al., 2013; Luu Chapman-Enskog expansion gives the incompressible Navier-stokes
et al., 2017; Ngoma et al., 2018; Tran et al., 2017). equations, and then the direct derivation of the macroscopic features of
V
the flow in the low Mach number limit, namely for Ma = max c
< < 1,
3.1.1. Solid phase modeling: Discrete Element Method (DEM) where Vmax is the maximum fluid velocity and c the lattice speed.
The DEM method, initially developed by Cundall and Strack (1979),

4
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 5. Scenario in the flooded (fully saturated) situation: dropout sinkhole and subsidence sinkhole, leading to inverted bowl and hourglass final shapes, respec-
tively.

from behavior at the grain scale to the erosion front kinetics, by going
through the mesoscale of force chains. In the present study, we focus on
cover-collapse sinkhole formation in a homogeneous soil layer (i.e.
without weaker zones). As presented in the previous section and illu-
strated in Fig. 5, we intend to reproduce the flood-induced collapse of a
cohesive sediment layer overlying karst caves or open fissures in the
bedrock, according to representative scenarios for the geological con-
text of the Orléans case study.
Our numerical model consists of a fully immersed cohesive granular
layer, at the bottom of which lies an orifice through which water flows
freely, assuming underneath a connection to the underground conduit
of a hydraulic network (Fig. 7). To induce erosion processes, we impose
Fig. 6. Illustration of the 2D coupled numerical method: DEM frictional contact a high enough pressure at the sample upper surface Pinlet, keeping the
model and LBM D2Q9 scheme (see text). outlet pressure at zero. For this study, we implement a cohesive as-
sembly of particles whose mean radius is r=1.5 mm. We vary the height
H and the length L of the 2D granular samples, corresponding to
3.1.3. Fluid-solid interaction
To couple the DEM and LBM methods, we assign a fluid or solid
status at each LBM grid point that is updated at each LBM time step
(with 2 DEM subcycles per LBM cycle). The space resolution for our
simulations is fixed so as to have a mean particle diameter equal to 10
LBM nodes. In the LBM code, when a velocity vector is directed towards
a node assigned as solid, we use a typical bounce-back algorithm for
which the fluid vector basically reverses its direction to obtain a no-slip
condition at any solid surface (grains and walls). To compute the hy-
drodynamic forces on grains, we refer to a relationship developed by
Bouzidi et al. (2001) based on a momentum-exchange between fluid
and solid nodes. Finally, the LBM calculation is implemented with a
hydraulic radius (dashed circles in Fig. 6) equal to 0.8 times the real
radius used in DEM. In this way, the fluid can flow between the disks of
the DEM assembly and we recover non-zero permeability in our 2D
granular configuration.

3.2. Application to sinkhole formation

Fig. 7. Numerical model for a sample with 7000 particles of mean radius
3.2.1. Numerical model setup r=1.5 mm (H=13.5 cm, L=40 cm, Pinlet=5 kPa and Δ=3 cm) at t = 0.6 s from
In a preliminary work, we reproduced the process of collapse in a the triggering of the pressure gradient. The particle cohesion number is
clogged underground conduit under high hydraulic gradients (Luu Coh=95, except in the dark grey zone of 1.5 cm thickness, where the cohesion
et al., 2017) and hence demonstrated the relevance of employing DEM- is 10 times greater in order to represent a more realistic boundary condition
LBM coupling to investigate the different scales involved in soil erosion, than a wall.

5
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Table 1 grains at the top. Subsequently, the destabilization front moves along a
Parameters used in the numerical simulations. direction opposite the flow direction (namely the “backward” direc-
Solid phase Fluid phase tion). We therefore observe a damage zone that gradually develops
upward, leading to the first scenario proposed in Fig. 5. The granular
−3
Mean radius, r 1.5⋅10 m Kinematic 4 10−5 m2/s medium self-organizes by creating force chains under compression.
viscosity, ν
Subjected to a traction mechanism, the erosion process takes place
Polydispersity, rmax/ 1.5 Lattice speed, c 10 m/s
rmin
through the progressive breaking of arches sustained by force chains.
Bond strength, C 10 N, 20 N Inlet pressure, [1–34] kPa Indeed, the total stress computation in Fig. 9 (right) gives a qualitative
Pinlet representation of how the force chains (displayed by the highest values)
Solid density, ρs 2.5⋅103 kg/m3 Fluid density, ρf 1 103 kg/m3 are involved in the erosion kinetics. The upward front advances in a
Geometry
series of vault-to-vault movements. Thus, our numerical model ap-
Length, L 23.7⋅10−2 m, 40⋅10−2
m proximately simulates the inverted bowl shape for the destabilization
Height, H [8.5–18.5]⋅10−2 m zone, hypothesized in Section 2. Note that the cavity formation could be
Orifice size, [15–75] ⋅10−3 m included in this backward erosion regime. It is actually likely that,
Conduit height 6⋅10−2 m
given a sufficient amount of simulation time, the cavity could expand
up to the surface by intermittent periods of grain evacuation.
Lastly, if the erosion parameters continue to increase (increasing
assemblies of 4835 to 9636 particles, while the size of the orifice Δ
Pinlet, increasing Δ and decreasing Coh), a continuous discharge regime
ranges from 10 to 50 times the mean particle radius. All the numerical
is recovered, similar to the subsidence sinkhole previously presented as
simulation parameters are given in Table 1 (note that the particle
the second speculative scenario in Figure5. In Fig. 10 (left), we observe
parameters are similar to those of usual 2D-DEM numerical studies).
that the evolution of the interparticle force network presents fracturing
To quantify soil cohesion, we define a dimensionless particle co-
initiated along two dislocation lines, forming a cone, and which pro-
hesion number as the ratio of the bond strength C to the particle's own
C pagates rapidly almost throughout the granular sample. The evacuation
buoyant weight, Coh = (ρ − ρ ) gS , where ρs − ρf is the submerged ap-
s f of grains through the orifice leads to the subsidence of the surface,
parent density, g is the gravitational acceleration and S is the particle's contrary to the previous backward erosion. This regime is very similar
surface (2D model). It is noteworthy that here particle cohesion consists to the hourglass and hopper, two classical granular matter problems.
of solid bonds with irreversible failure, and is different from the mac- But in our case, we added interstitial fluid flow and interparticle co-
roscopic cohesion classically defined through the Mohr-Coulomb failure hesion. This regime of downward migration of cohesive particles
criterion. The range of Coh studied is equivalent to a weakly cemented (which is not usual in geotechnics, see previous Section in 2.3) has been
granular matter. Regarding flow regime, the Reynolds number can be recently studied numerically using CFD-DEM coupling (Liu et al.,
measured within the bottom conduit before the detachment of the 2018), but no experimental work has been dedicated to this phenom-
VΔ
grains, such that Re = ν , where V is the fluid mean velocity and ν the enon to our knowledge.
fluid kinematic viscosity. Globally, in the granular assembly and out-
side, V is of the same order of magnitude as the maximum fluid velocity 3.2.3. Phase diagrams
Vmax, ranging from 0.01 to 1 m/s. For reasons of computation time, the Figs. 11 and 12 present phase diagrams that locate the domains of
particle size is taken substantially larger than in reality. Therefore, we existence of cavity formation, dropout sinkhole and subsidence sink-
use a higher viscosity than water to simulate realistic flow regimes by hole, according to inlet pressure Pinlet, orifice size Δ, sample height of H,
obtaining similar Reynolds numbers to those in the field. Thus, all the sample length L and cohesion number Coh. Fig. 11 shows two diagrams
following simulations reproduce flow regimes that are laminar and displaying Pinlet as a function of Δ for a fixed Coh, a fixed L, and two
incompressible, with an Re from 30 to 600 and a Mach number different sample heights of H=13.5 cm (a) and H=8.5 cm (b). For both
V
Ma = max c
smaller than 0.1. cases, the erosion process successively jumps from cavity to dropout
sinkhole and ultimately to subsidence sinkhole when Pinlet and Δ are
3.2.2. Phenomenology increased. It can be seen that, for the range studied, the location of the
We performed a parametric study using our DEM-LBM code to in- frontier between two regimes depends on the height H of the granular
vestigate erosion processes according to two different micromechanical sample. In an attempt to rationalize this dependency, we plotted (not
viewpoints: from the bonds within the cohesive granular sample and shown here) the dimensionless pressure, dividing Pinlet by the hydro-
from the local stress field. The following figures (Figs. 8, 9, 10) present static pressure ρfgH, as a function of different dimensionless orifice sizes
typical simulations of microstructure evolution over time regarding: on (using r, H or L). But in each case the regime transitions kept being
the left side, the interparticle normal force network composed of po- influenced by H. In Fig. 12a, we nevertheless found empirically that
sitive (compressive) and negative (tensile) normal forces. On the right multiplying Pinlet by H leads to an approximate grouping of the different
side, the total stress tensor for each grain, as the contribution of both erosion regime frontiers, regardless of the sample geometry (for various
the stress exerted by the neighboring particles in contact and the hy- H and L) and for a given cohesion state. Although this result could have
draulic stress of the surrounding fluid. We choose the stress tensor practical implications, we do not yet have any physical explanation for
determinant to quantify this magnitude. it. Finally, Fig. 12b indicates that for a doubled Coh, a higher Pinlet and
The first erosion regime is illustrated in Fig. 8. From the bottom Δ are necessary to pass from dropout to subsidence sinkholes, and there
orifice we observe an evacuation of grains that stops early, leaving a is no longer any cavity formation within the same range of parameters.
stable cavity. This is a marginal regime also simulated in our previous
work on clogged underground conduit collapse [21]. This vault cavity 3.2.4. Spatio-temporal diagrams
recalls quasi-static intermittent void formations in quasi-static and non- In this numerical study, we developed a post-processing tool to
flooded conditions, commonly described in geotechnical works carry out a spatio-temporal analysis on the basis of our simulations.
(Sowers, 1996). Fig. 13 presents a typical image of the interparticle force network in its
Then, when increasing the hydraulic load (Pinlet), or enlarging the initial state. Using ImageJ, an open-source software application, we can
bottom orifice (Δ), or decreasing the interparticle cohesion (Coh), we select the main eroded zone corresponding to the orifice width (on the
observe the expected backward erosion implied in the dropout sinkhole. left image) to plot averaged spatio-temporal diagrams (on the right).
Fig. 9 (left) shows how interparticle bonds progressively break from the Thus, the vertical axis of a diagram displays the mean pixel density
orifice to the surface. In the early stages, there is no disturbance of the calculated in the selected zone and the horizontal axis is time. Using

6
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 8. Cavity formation for a sample with 7000 particles of mean radius r=1.5 mm (H=13.5 cm, L=40 cm, Pinlet=0.5 kPa, Δ=4.5 cm, and Coh=95) with evolution
of the microstructure over time: normal interparticle force (left) and total stress determinant (right).

this representation, clearly shows that the two erosion regimes mainly loading. To this end, we perform a systematic procedure to characterize
differ by their front propagation directions. The spatio-temporal dia- the maximum width of the eroded zone. On the spatio-temporal dia-
gram related to the dropout sinkhole shows an almost constant surface gram presented previously, we detect a characteristic time t0 at which
level with time, while the backward erosion advances from the bottom. the surface starts to cave-in, as indicated in Fig. 14a. At that moment
Conversely, the subsidence sinkhole displays an almost linear decrease t = t0, we apply the Plot Profile function of ImageJ that displays a two-
of the surface level with time. This difference in the induced ground dimensional graph of the intensities of pixels along a line. The x-axis
motion (observed from the surface) is indeed that commonly agreed by represents the horizontal distance and the y-axis the vertically averaged
the natural hazards community to distinguish dropout from subsidence pixel intensity (Fig. 14b). Thus, we can estimate the so-called collapse
sinkholes (Sowers, 1996). width λ with acceptable measurement uncertainty.
By applying this image processing for typical simulations of the
dropout sinkhole regime, we obtain the graph shown in Fig. 15 (left).
3.3. Predictive analysis for dropout sinkholes The collapse width λ is plotted as a function of the height of the
granular sample H. The main finding is that λ is roughly equal to the
The objective of our research is to improve the prediction of where, cohesive soil thickness H (slope equal to 1), accounting for a moderate
why and when a sinkhole will occur. From the risk assessment view- discrepancy and irrespectively of H, Pinlet, Δ and Coh. When considering
point, the dropout sinkhole is the most critical one because cavity second-order trends, we remark that the higher Pinlet and Δ are, the
widening cannot be viewed directly from the surface. A relevant aspect higher λ becomes. When Coh is doubled, λ is either slightly smaller or
of the micromechanical modeling we propose is precisely to give access equal within the error bar.
to the internal evolution of the process. Therefore, we present here a With a first-order approximation, the simulations of the dynamic
qualitative first-order predictive analysis based on the numerical si- process in a flood situation therefore predict a linear relationship be-
mulations of the backward erosion scenario. tween the collapse size λ and the cover thickness H, according to
As mentioned in the previous section dedicated to the field sce- Sower's prediction but with a higher slope of 1 instead of 2/3. In ad-
narios, Sowers proposed the geotechnical approach of considering a dition, we find that this result does not depend on the other parameters
potential collapse diameter of 2/3 the height of the overlying soil layer (Pinlet, Δ and Coh), therefore providing an interesting relationship for
(Sowers, 1996). This approach concerns the quasi-static growth of a applications. In contrast, as expected, Fig. 15 (right) indicates that the
void in a cohesive but non-submerged soil. Using our simulations, we erosion kinetics depends on all the parameters studied. Indeed, by
propose to estimate the extension of bond breakage in a flood situation, plotting t0 the time at which the surface starts becoming destabilized
when backward erosion occurs in a saturated soil under hydraulic

7
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 9. Dropout sinkhole for a sample with 7000 particles of mean radius r=1.5 mm (H=13.5 cm, L=40 cm, Pinlet=1.0 kPa, Δ=4.5 cm, and Coh=95) with evolution
of the microstructure over time: normal interparticle force (left) and total stress determinant (right).

(defined in Fig. 14a) as a function of H, we observe a much broader together with liquid paraffin. The amount of paraffin is chosen low
spread of data than for λ. For a fixed H, whereas the maximum dis- enough to prevent the complete filling of voids between particles, in the
crepancy for λ corresponds to a doubled value, the scatter of t0 can sense that the liquid far from saturates the granular sample. We then
reach up to a factor 6. A more systematic study would be necessary to paid particularly attention to the homogeneous mixing of the sand with
capture the dependency of t0 on the different parameters. Qualitatively, hot liquid paraffin. When the mixture is cooled and the paraffin has
we observe that the higher Pinlet and Δ is, the faster the erosion. In- solidified, we obtain an artificial cohesive granular material that can be
versely, the greater the H, the slower the erosion kinetics. In these si- considered as a weakly cemented sand.
mulations, particle cohesion seems to have little influence. Note that The experimental setup is pictured in Fig. 16. The cohesive sand is
the estimation of the characteristic collapse time t0 here gives only prepared in a specific sample cell open at the top and open only by a
qualitative but not quantitative information on time scales. The next circular orifice at the bottom. For this study, we tested an orifice dia-
step would be to calibrate using real soil parameters to obtain a more meter between 1.5 cm and 2.5 cm. By filling this cell, we produced soil
realistic prediction of sinkhole dynamics, which would greatly improve samples whose height varied between 10 cm and 12 cm, with a fixed
risk management. Linking the numerical modeling illustrated here to length of 20 cm and a fixed width of 5.5 cm. This sample cell was then
geotechnical assessments will be a significant step towards character- introduced into a larger (external) cell connected to a gear pump. This
izing sinkhole hazards. closed-loop system induced a controlled flow rate to continuously erode
the granular sample with a downward water flow. During the erosion
4. Experimental modeling experiment, the sand was assumed to be evacuated through the output
conduit connected to the upper reservoir. In all the experiments per-
In addition, we performed an exploratory experimental study to formed, we could see through the transparent plastic tube that eroded
enrich the previous numerical results. An original methodology has sand flowed out of the orifice without clogging the pipe. A primary
recently been developed to produce artificial cohesive granular media saturation phase is crucial before erosion is triggered. Contrary to the
similar to our numerical granular assemblies (Brunier-Coulin, 2016). flow direction for erosion (“input” to “output”), we allowed the water
to rise very slowly from the bottom with the air being removed up-
4.1. Materials and setup wards. In this way, we were able to immerse all the pores in the
granular matter.
For the present study, we use Hostun sand HN31 (mean diameter of
240 μm), which is a reference sand used in many experimental geo-
technical studies. Solid bridges are created by “sticking” sand grains

8
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 10. Subsidence sinkhole for a sample with 7000 particles of mean radius r=1.5 mm (H=13.5 cm, L=40 cm, Pinlet=9.6 kPa, Δ=4.5 cm, and Coh=377) with
evolution of the microstructure over time: normal interparticle force (left) and total stress determinant (right).

4.2. Results paraffin concentrations, i.e. around 0.3% per weight for both examples.
For higher values, the maximum pump flow rate was not high enough
As can be seen in Fig. 17, this rather simple setup nevertheless al- to erode the soil. For lower concentrations of paraffin, the cohesive
lowed us to simulate two regimes very similar to those observed in our sample was too easily weakened during its preparation and introduc-
numerical study and consistent with the field scenarios. However, re- tion in the setup, leading to fractures and even collapse when the
garding sand cohesion, we were unable to explore a wide range of eroding flow was imposed. Moreover, both processes were obtained

(a) (b) 8000

8000
Subsidence
7000 Subsidence 7000

6000 6000

5000 5000
Pinlet (Pa)
Pinlet (Pa)

4000 4000
Dropout
3000 3000
Dropout 2000
2000

1000 1000
Cavity Cavity
0 0
0 0.03 0.06 0.09 0.12 0 0.03 0.06 0.09 0.12
(m) (m)

Fig. 11. Phase diagrams related to the inlet pressure Pinlet and the orifice size Δ, (a) for Coh=95, H=13.5 cm and L=40 cm. (b) Idem for Coh=95, H=8.5 cm and
L=40 cm.

9
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

(a) 1200 (b) 1200

Subsidence
1000 1000 Subsidence

Pinlet H (Pa.m)
800 800
Pinlet H (Pa.m)

600 600
Dropout

400 400
Dropout

200 Cavity 200

0 0
0 0.03 0.06 0.09 0.12 0 0.03 0.06 0.09 0.12
(m) (m)

Fig. 12. Phase diagrams related to the inlet pressure Pinlet multiply to H and the orifice size Δ, (a) for Coh=95, for different sample heights H=8.5 cm (square),
13.5 cm (circle) and 18.5 cm (triangle) when L=40 cm, and H=12.8 cm, L=23.7 cm (diamond). (b) Idem for Coh=191, L=23.7 cm and H=13.5 cm (circle).

with the same flow rate. Inhomogeneity in the sample was also iden-
tified. This could explain why we observed large aggregates and falling
blocks. Therefore, as currently framed, the experimental conditions
allowed us neither to precisely control the parameters nor to select one
of the two regimes. Improvement of the protocol is in progress.
In spite of this, we were able to obtain an interesting result from the
dropout sinkhole process. As indicated in Fig. 17, we estimated a
characteristic damage size λ. As can be seen in Fig. 15 (red losange),
this measurement is in good agreement with the prediction of Sowers
(lower than predicted by our simulations, as discussed in 3.3), giving
relevancy to this exploratory experimental study. In future work, to
ensure better comparability with our numerical simulations, certain
tests could be performed by varying the sand diameter to the cell width
ratio in order to study the 3D effects. The next step will also focus on
the calibration of parameters such as artificial interparticle cohesion
using purely mechanical tests (traction test, direct shear), in order to
make more quantitative comparisons between the experiments and si-
mulations.

5. Conclusion and perspectives

Sinkhole occurrences represent a major geohazard that remains

poorly assessed, especially in exceptional flood situations. The me-
teorological event occurring in the Orléans area in spring 2016, is a Fig. 14. Dropout sinkhole case: (a) snapshot of interparticle bonds at the
relevant example of a huge increase in the number of cavity collapses. characteristic time t0 indicated on the averaged spatio-temporal diagram. (b)
This paper presented a physical model of hydro-mechanical instabilities Plot Profile (using ImageJ software) at t0 and collapse width λ estimation (see
under hydraulic forcing on a cohesive soil layer underlain by a karstic text).
cavity. The numerical study involving a coupled LBM-DEM method
provided insight into the micromechanism of particle bond breakage

Fig. 13. Typical spatio-temporal diagrams averaged within the central orifice zone (red frame), and showing typical evolutions for both dropout and subsidence
sinkholes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

10
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 15. (a) Collapse width λ vs granular sample height H, for various inlet pressures Pinlet, orifice size Δ and particle cohesion Coh. Linear fit giving a slope equal to 1.
The grey area indicates the 95% confidence bounds obtained on the slope. (b) The corresponding characteristic collapse time t0.

during internal erosion processes. By varying hydraulic pressure, un- erosion process, consistent with classical geotechnics estimations of
derground conduit size and the particle cohesion of the granular media, hypothetical sinkhole size. Finally, the first results of an experimental
we identified two different erosion kinetics (excluding cavity forma- model demonstrated its capacity to produce similar sinkholes, as pre-
tion), namely upward cavity development and downward granular dicted by the field scenario and numerical simulations.
discharge, leading to the formation of dropout and subsidence sinkholes Following these results, one of the immediate future goals will be to
(in cohesive soil), respectively. On the basis of a parametric study, we integrate calibration in our numerical model. Further comparison be-
plotted a phase diagram to delimit the respective domain of each re- tween both numerical-experimental approaches will include tran-
gime. Then, we performed systematic measurements, directly from our scribing the cohesion achieved in the laboratory tests (typically traction
simulations, of the eroded zone when a dropout sinkhole occurred. We and shear tests) to the particle bond strength implemented in the co-
found a linear relationship between the characteristic collapse width hesion model of our code. Linking this multi-scale numerical model to
and the height of the cohesive soil overlaying the underground conduit. geotechnical in situ conditions will be a significant step in character-
This result demonstrated the ability of our calculations to reproduce an izing sinkhole hazards and risks.

Fig. 16. Experimental setup. Left: Sample cell. Right: hydraulic closed-loop (see text).

11
L.-H. Luu, et al. Engineering Geology 260 (2019) 105249

Fig. 17. Typical experimental pictures of the two different erosion regimes. For the dropout and subsidence sinkhole, the sample height is equal to 10 cm and 12 cm,
and the orifice diameter around 1.5 cm and 2.5 cm, respectively. The estimation of the maximum damage λ in the dropout sinkhole case is indicated.

Acknowledgements Lominé, F., Scholtès, L., Sibille, L., Poullain, P., 2013. Modeling of fluid-solid interaction
in granular media with coupled lattice Boltzmann/discrete element methods:
Application to piping erosion. Int. J. Numer. Anal. Methods Geomech. 37, 577–596.
We thank André Benjamin (internship student) for his contribution Lorain, J.M., 1973. La géologie du calcaire de Beauce. Le Calcaire de Beauce. Bulletin de
to the experimental study. This research was funded by PERCIVAL liaison des laboratoires des Ponts et Chaussées. 4.
(PErception des Risques effondrements liés aux Cavités associés aux Luu, L.-H., Philippe, P., Noury, G., Perrin, J., Brivois, O., 2017. Erosion of cohesive soil
layers above underground conduits. In: EPJ Web of Conferences. 140 EDP Sciences,
Inondations en VAL de Loire), French Regional Research Program 09038.
(2018-2020) attributed by the Centre-Val de Loire region. Moreau, J., 2002. Les gouffres du nord d’Orléans. Groupe spéléologique orléanais.
Bulletin 7, 1–28.
Ngoma, J., Philippe, P., Bonelli, S., Radjaï, F., Delenne, J.Y., 2018. Two-dimensional
References numerical simulation of chimney fluidization in a granular medium using a combi-
nation of discrete element and lattice Boltzmann methods. Phys. Rev. E 97 (5),
Baryakh, A.A., Fedoseev, A.K., 2011. Sinkhole formation mechanism. J. Min. Sci. 47 (4), 052902.
404–412. Noury, G., Perrin, J., Luu, L.H., Philippe, P., Gourdier, S., 2018. Role of floods on sink-
Beck, B., 2012. Soil piping and sinkhole failures. In: Encyclopedia of Caves (Second holes occurence in covered karst terrains: Case study of Orléans area (France) during
Edition) pp. pp. 718–723. the 2016 meteorological events and perspectives for other karst environments. In:
Benahmed, N., Canou, J., Dupla, J.C., 2004. Initial structure and static liquefaction 15th Multidisciplinary Conference on Sinkholes and the Engineering and
properties of sand. Comptes Rendus Mecanique 332 (11), 887–894. Environmental Impacts of Karst.
Bonelli, S., 2013. Erosion in Geomechanics Applied to Dams and Levees. John Wiley Sons. Perez, A.L., Nam, B., Chopra, M., Sallam, A., 2017. Understanding Florida’s Sinkhole
Bouzidi, M.H., Firdaouss, M., Lallemand, P., 2001. Momentum transfer of a Boltzmann- Hazards: Hydrogeological Laboratory Study. In: Geotechnical Frontiers 2017, pp.
lattice fluid with boundaries. Phys. Fluids 13 (11), 3452–3459. 508–518.
Brinkmann, R., Parise, M., Dye, D., 2008. Sinkhole distribution in a rapidly developing Perrin, J., Cartannaz, C., Noury, G., Vanoudheusden, E., 2015. A multicriteria approach to
urban environment: Hillsborough County, Tampa Bay area, Florida. Eng. Geol. 99 karst subsidence hazard mapping supported by weights-of-evidence analysis. Eng.
(3–4), 169–184. Geol. 197, 296–305.
Brunier-Coulin, F., 2016. Etude des mécanismes élémentaires de l’érosion d’un sol cohésif. Radjaï, F., Dubois, F., 2011. Discrete-Element Modeling of Granular Materials. Wiley-Iste,
(Doctoral dissertation, Aix-Marseille). pp. 425.
Cuellar, P., Philippe, P., Bonelli, S., Benahmed, N., Brunier-Coulin, F., Ngoma, J., Radjai, Rawal, K., Hu, L.B., Wang, Z.M., 2017. Numerical investigation of the geomechanics of
F., 2015. Micromechanical analysis of the surface erosion of a cohesive soil by means sinkhole formation and subsidence. In: Geotechnical Frontiers 2017, pp. 480–487.
of a coupled LBM-DEM model. In: Proceedings of the International Conference on Sellmeijer, H., de la Cruz, J.L., van Beek, V.M., Knoeff, H., 2011. Fine-tuning of the
Particles, Barcelona, Spain, pp. 28–30. backward erosion piping model through small-scale, medium-scale and IJkdijk ex-
Cundall, P.A., Strack, O.D.L., 1979. The development of constitutive laws for soil using periments. Eur. J. Environ. Civ. Eng. 15 (8), 1139–1154.
the distinct element method. Numerical Methods Geomech. 1, 289–317. Shalev, E., Lyakhovsky, V., 2012. Viscoelastic damage modeling of sinkhole formation. J.
De Waele, J., Gutiérrez, F., Parise, M., Plan, L., 2011. Geomorphology and natural hazards Struct. Geol. 42, 163–170.
in karst areas: a review. Geomorphology 134 (1–2), 1–8. Silvani, C., Désoyer, T., Bonelli, S., 2009. Discrete modelling of time-dependent rockfill
Delenne, J.Y., El Youssoufi, M.S., Cherblanc, F., Bénet, J.C., 2004. Mechanical behaviour behaviour. Int. J. Numer. Anal. Methods Geomech. 33 (5), 665–685.
and failure of cohesive granular materials. Int. J. Numer. Anal. Methods Geomech. 28 Sowers, G.F., 1996. Building on Sinkholes: Design and Construction of Foundations in
(15), 1577–1594. Karst Terrain. American Society of Civil Engineers.
Ford, D., Williams, P.D., 2013. Karst Hydrogeology and Geomorphology. John Wiley Talon, L., Bauer, D., Gland, N., Youssef, S., Auradou, H., Ginzburg, I., 2012. Assessment of
Sons. the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous
Gutiérrez, F., Parise, M., De Waele, J., Jourde, H., 2014. A review on natural and human- media. Water Resour. Res. 48 (4).
induced geohazards and impacts in karst. Earth Sci. Rev. 138, 61–88. Tao, X., Ye, M., Wang, X., Wang, D., Pacheco Castro, R., Zhao, J., 2015. Experimental and
Hyatt, J.A., Jacobs, P.M., 1996. Distribution and morphology of sinkholes triggered by Numerical Investigation of Sinkhole Development and Collapse in Central Florida.
flooding following Tropical Storm Alberto at Albany, Georgia, USA. Geomorphology Tran, D.K., Prime, N., Froiio, F., Callari, C., Vincens, E., 2017. Numerical modelling of
17 (4), 305–316. backward front propagation in piping erosion by DEM-LBM coupling. Eur. J. Environ.
Iverson, R.M., 2000. Landslide triggering by rain infiltration. Water Resour. Res. 36 (7), Civ. Eng. 21 (7–8), 960–987.
1897–1910. Waltham, T., 2008. Sinkhole hazard case histories in karst terrains. Q. J. Eng. Geol.
Kuniansky, E.L., Weary, D.J., Kaufmann, J.E., May 2016. The current status of mapping Hydrogeol. 41 (3), 291–300.
karst areas and availability of public sinkhole-risk resources in karst terrains of the Wilson, G.V., Rigby, J.R., Dabney, S.M., 2015. Soil pipe collapses in a loess pasture of
United States. Hydrogeol. J. 24 (3), 613–624. Goodwin Creek watershed, Mississippi: role of soil properties and past land use. Earth
Lallemand, P., Luo, L.S., 2000. Theory of the lattice Boltzmann method: Dispersion, Surf. Process. Landf. 40 (11), 1448–1463.
dissipation, isotropy, Galilean invariance, and stability. Phys. Rev. E 61 (6), 6546. Yu, D., Mei, R., Luo, L.S., Shyy, W., 2003. Viscous flow computations with the method of
Liu, P., LaMarche, C.Q., Kellogg, K.M., Hrenya, C.M., 2018. A square force cohesion lattice Boltzmann equation. Prog. Aerosp. Sci. 39 (5), 329–367.
model and its extraction from bulk measurements. AICHE J. 64 (7), 2329–2339.

12