NUMERICAL SIMULATION OF KARST SOIL CAVE EVOLUTION

arching effect in the overlying soil layer. By analyzing the plastic zone of the soil layer, it was found that, in rigid clay, arch roof collapse and tensile failure are the major events that lead eventually to the barrel-shaped or bottle-shaped forms of collapsed pits. In loose soil, shear failure of the arch toe is the major event that eventually leads to the taper-shaped or bowl-shaped form of collapsed pits. Generally speaking, stability and size of soil caves can be determined using the three variables of low shear stress area, equal settlement plane, and plastic zone discussed above. The numerical simulation of this study is valuable to the monitoring and assessment of sinkhole occurrence.


Introduction
A cover-collapse sinkhole, a roughly circular hole in the ground that is variable in depth, is one of the most adverse geological phenomena in an area underlain by limestone bedrock.The collapse feature is a result of the movement of soil or other related materials carried by water down into voids either in the limestone bedrock or within the soil profile.This study aims to understand the formation and development of voids in the soil profile (the karst soil cave) related to the occurrence of cover-collapse sinkholes, which has been a long-lasting challenge for various reasons, including the hidden development of the soil caves and complexity in karst hydrogeology and geological conditions.The main approaches of studying sinkhole soil caves are as follows:

Abstract
This study is focused on numerical simulation of the formation and development of karst soil caves related to cover-collapse sinkholes.The so-called 'karst soil cave' refers to the caves formed in the soil layers above bedrock of sinkhole regions.Because the soil caves are formed and developed under groundwater seepage, studying groundwater level changes can help understand soil cave development and collapse.Based on the improved Terzaghi loosening pressure theory and using excess pore water pressure, two kinds of critical groundwater level decline are discussed.The first, denoted as ΔH 0 , is the critical groundwater level decline related to soil cave formation and evolution; the other one, denoted as ΔH T , is the critical groundwater decline related to cave roof collapse.After a soil cave is formed, its evolution can cause uneven displacement and stress redistribution in the overlying soil layer.The process of soil cave expansion can be understood by investigating the change in displacement and stress.Numerical simulation of the vertical displacement using FLAC3D shows that the maximum vertical displacement occurs at the arch roof of the soil cave and that the displacement can cause tensile failure of the arch roof.The simulated soil layer displacement is used to determine the soil depth disturbed by the cave by delineating the planes of equal settlement.Analyzing the simulated shear stress shows that the maximum shear stress occurs at the arch toes and causes shear failure.On the other hand, the zones of low shear stress can be used to evaluate existence of

Trigger Factor Method
In China, groundwater level decline is the major triggering factor of sinkhole collapse, and the decline plays a key role in soil cave formation and evolution (Figure 1).Therefore, in areas prone to sinkhole collapse, monitoring changes in groundwater levels can be used for early warning of sinkhole collapse.When water level change exceeds a critical value, sinkhole formation is expected.To implement this requires determining the critical conditions of karst collapse (Lei 2010, Yan 2014), and this is one of the motivations for this study.
Another motivation for this study is to understand sinkhole evolution caused by uneven displacement and stress redistribution of the soil overlying the sinkhole.The displacement and stress variation is simulated numerically in this study using FLAC3D, a threedimensional finite difference program.Figure 2 is an example computational grid of FLAC3D used to simulate sinkhole evolution.

Geophysical Methods
Ground-penetrating radar (GPR) is a relatively popular geophysical method for investigating sinkholes.The GPR relies on the reflection of a high frequency (25-1,000 MHz) electromagnetic (EM) pulse from layer contacts and other "anomalies" (boulders, cavities, utilities, tanks, etc.).GPR has become an important non-invasive technique for rapid geophysical mapping of the shallow subsurface.(Benson, 1987;Witten and Calvert, 1999;Martin-Crespo and Gomez-Ortiz, 2007).However, because of the signal attenuation, depth of penetration is limited.For example, subsurface clay always attenuates GPR signals.This occurs often because most cover-collapse sinkholes are related to clayey overburden on cavity-laden limestone bedrock.Other geophysical methods are available, such as Direct current (DC) resistivity techniques (Panno et al., 1994;Batayneh and Al-Zoubi, 2000), micro-gravity (Butler, 1984) and electromagnetic method (Kaspar and Pecen, 1975).Although these techniques have been widely used in karst regions, they cannot be used for effective monitoring and early warning of sinkhole collapse.

Remote Sensing Methods
Remote sensing methods, such as InSAR (Interferometric Synthetic Aperture Radar), provide satellite images of the Earth's surface that can be combined to show subtle movements, i.e., deformation, of the ground surface.The ground deformation can be used to monitor and provide early warning of sinkhole collapse.However, remote sensing methods are unsuitable to real-time monitoring (Paine et al., 2009).Similarly, LiDAR (Light Detection and Ranging) technology is unsuitable to real-time monitoring, although it provides high density data which can be processed to produce a high resolution topographic map (Doctor and Young, 2013;Shaw-Faulkner et al., 2013).

Photoelectric Monitoring Methods
These methods include optical fiber sensing (BOTDR and OTDR) and coaxial cable sensing (TDR).By burying fiber or coaxial cable in soil, soil deformation and damage creates strain on the cables, and the strain can be measured by special instruments.This makes it possible for real-time monitoring.Photoelectric sensors have played the role of precise monitoring and early warning for small-scale engineering sites.However, photoelectric sensors have not been applied to a largescale site (Dowding 2003, Guan et al. 2013).where D is hole diameter, K is the lateral earth pressure coefficient, Z is the depth of soil and the roof panel, and is the angle of internal friction.Equation (3) shows that, when the soil hole diameter is small, the σ z is negative.Under this condition, there exists a collapse resistance, and the soil cave is temporarily stable.Analysis of the above equations reveals the following: 1.When the overlying layer of the soil cave is sand and the cohesion c is very small, ΔH 0 is small.As a result, small fluctuation of groundwater level may cause formation of soil cave.In addition, σ z in a small soil cave is already positive, and it leads to emergence of roof collapse due to raveling of soil into the void space in the underlying bedrock 2. When the overlying layer is clay and cohesion c is big, soil cave can be in the temporarily stable state.When groundwater level declines more than ΔH 0 , soil cave starts to grow upward to the ground surface until a new stable state is formed or collapse occurs.
3. When groundwater level dropped more than ΔH T , soil cave roof is unstable, and collapse occurs.When the soil hole diameter is small, ΔH T > ΔH 0 when the soil hole become larger, ΔH T becomes smaller.
For all the cases of sinkhole formation and collapse, the overlying layer is subject to displacement and stress changes.It is therefore necessary to understand the changes.FLAC3D is used in this study to simulate soil displacement and stress in the stable stage of sinkhole development.The FLAC3D-based model uses a three-dimensional domain with the dimension of a×b×c, and the dimension can be adjusted to meet specific project needs.As shown in Figure 2, the model grid consists of tetrahedral blocks.The bottom boundary (the soil rock interface) of the domain is set as a fixed constraint boundary, and the other boundaries are set as the one-way boundary.The soil cave roof is always an arch, and generalized in this study as a hemisphere with the height of H and the roof span of D. The size parameters can be adjusted for different projects by using the FISH language of FLAC3D programming.The FLAC3D modeling of displacement and stress uses the Mohr-Coulomb failure criterion, and the associated parameters are listed in Table 1.

Soil Cave Development Process and Numerical Model
The soil cave formation and development is closely related to the properties and structures of soil, water, and rock.Figure 1 illustrates the process of sinkhole formation and development.First, the soil cave formation requires a karst cave opening upward and having sufficient space to accommodate falling soil.Under gravity and groundwater seepage forces, the soil layer breaks and enters the karst cave, and this is the formation of soil cave.Afterward, the soil cave can quickly expand upward.When shear stress is less than cohesive strength of the soil, the soil cave stops developing and remains in a temporarily stable state.When the dynamic conditions change, the soil sinkhole will continue to expand until the collapse of the soil cave roof, as shown in Figure 1.where c is soil cohesion, g is the gravitational acceleration, ρ w is the water density, and σ t is soil tensile strength.
For the collapse of a soil cave, when the summation of excessive pore water pressure and vertical pressure is positive, i.e., σ w + σ z > 0, the soil cave roof collapses, and the critical water level drop for cave roof collapse is 2 The soil sinkhole vertical pressure σ z without overlying load can be derived using the Terzaghi loose pressure theory improved by Kezdi (1975) for a circular cavern.
It is expressed as The simulation procedure is as follows: 1. Conduct initial stress analysis for the entire modeling domain before the soil cave is formed.
2. Generate the cave (with the height of H and span of D) using the nill command of the FLACD FISH language.
3. Conduct stress analysis to simulate stress and strain of the new model domain with the cave.
4. Use the nill command to expand the cave size and simulate the sinkhole evolution and development.

Analysis of Vertical Displacement
In order to analyze the displacement changes during the soil cave expansion and evolution, the vertical displacement is simulated using FLAC3D and the contour lines are shown in Figure 3.The simulation results indicate that the soil settlement is uneven due to the emergence of the soil cave; the maximum vertical displacement always occurs in vault.The vertical displacement can cause tensile failure, extending from the vault to both sides.Correspondingly, ground subsidence also occurs, and the largest subsidence generally appears above the vault.
The vertical displacement of the overlying layer can be used to determine the position of the "equal settlement plane".The plane is defined as the position where uneven settlement begins.In China, the value of allowable settlement is 30mm, and is used in the community of geological engineering.This value is used to determine the "equal settlement plane", i.e., lines m-n marked in Figure 3.With the soil cave expanding upward, the  position of "equal settlement plane" rises accordingly.
When the plane reaches the ground, the whole overlying layer subsides, and the soil cave is in the critically unstable stage for sinkhole collapse to occur.

Analysis of Shear Stress
Using the FLAC3D FISH language, the contours of maximum shear stress are plotted in Figure 4.The maximum shear stress is half of the difference of the maximum principal stress and the minimum principal stress, i.e., τ max =σ 1 −σ 3 /2.Figure 4 shows the phenomenon of stress concentration at the arch toes.If the maximum shear stress near the toes reaches the soil strength, then the soil cave is subject to shear failure, and the damage spreads from the toes to the arch vault.
During the simulation process, it is found that a "low value zone" of shear stress appears at the top of the soil cave, where the shear stress contours are depressed downward on top of arch (Figure 4a, b).This phenomenon can be explained by the soil arch effect (Terzaghi, 1943;Handy, 1985).Due to the uneven displacement of the overlying layer, the stress of the entire domain is redistributed, and the stress applied to the arch is spread out to the toes and the surrounding media.
As a result, the maximum principal stress, σ 1 , decreases, and the minimum principal stress, σ 3 , increases, causing diminution of the maximum shear stress, τ max .Therefore, the shear stress "low value zone" can be used to determine the stabilization of the soil cave.The "low value zone" also evolves with the cave development.As shown in Figure 4, as the soil cave extends upward, the "low value zone" becomes gradually smoother (Figures 4a and 4b), and eventually disappears when the soil cave develops to a certain extent (Figure 4c).At this moment, the soil arch effect disappears, and the soil cave is in the unstable state and sinkhole collapse can occur.

Simulation of Plastic Zone
FLAC3D can be used to simulate the plastic zone for examining the area of potential failure regions around the soil cave.In the rigid clay layer, the plastic zone is mainly above the arch of the soil cave, and becomes bigger with the cave development (Figure 5, subplots a1, b1 and c1).Under this circumstance, the overlying soil layer is subject to tensile failure, and the vault collapses to form a barrel-shaped or bottle-shaped pit.At the same time, due to stress release and the existence of the free face, the soil near the upper arch toes inclines and slides to the pit center.Concentric tension cracks also form on the ground toward the collapse pit.
In the loose sand layer, the plastic zone appears near the arch toes, and develops along the tangential direction of the arch as the soil cave further develops (subplots a2, b2, and c2 of Figure 5).Therefore, the cave collapse starts from the shear failure at the arch toes.Under this circumstance, the entire overlying layer is broken, and taper-shaped or plate-shaped collapse pit is formed.In the early stage of this process, the plastic zone is formed only in the periphery area of the soil cave.With the development of the soil cave, the plastic zone occurs throughout the entire overlying layer.When the layer falls within the plastic zone, ground collapse occurs.Therefore, the stage of soil cave formation can be determined by examining the changes of the plastic zone.

Figure 1 .
Figure 1.Sketch of the process of soil cave development.

Figure 2 .
Figure 2. Model domain and grid of soil cave simulation using FLAC3D.
Groundwater, as the most active factor of sinkhole collapse, plays a key role in the formation and evolution process of karst soil caves.Based on soil sinkhole experiments,Jiang (1998) proposed that, when the change of karst cavity pressure reaches a certain value, soil failure occurs.Since karst cavity pressure and karst aquifer water pressure are subject to the same changes, formation of soil cave and collapse can be studied by investigating karst aquifer water level changes.Based on the relations of excessive water pressure, σ w ≈ρ w gΔH, and σ w ≥[σ t ]≈c/3,Wan (2003) developed the critical water level decline for soil sinkhole cave formation and continuous development 1

Figure 3 .
Figure 3. Diagram of FLAC3D-simulated vertical displacement with the cave dimension of (a) 5m, (b) 5m, and (c) 10m.D is the cave span.The overlying soil is clay.The lines of "m-n" denote the "equal settlement plane".

Figure 4 .
Figure 4. Distribution of FLAC3D-simulated shear stress with the cave dimension of (a) 5m, (b) 5m, and (c) 10m.D is the cave span.The overlying soil is clay.