**Pre-conference publication**
22 Apr 2020

**Pre-conference publication** | 22 Apr 2020

# A New Software to Model Earth Fissure Caused by Extensive Aquifer Exploitation and its Application to the Guangming Village Case, China

Yueting Li Matteo Frigo Yan Zhang Lin Zhu Massimiliano Ferronato Carlo Janna Xulong Gong Jun Yu Pietro Teatini and Shujun Ye

^{1},

^{2},

^{3},

^{4},

^{2},

^{2},

^{3},

^{3},

^{2},

^{1}

**Yueting Li et al.**Yueting Li Matteo Frigo Yan Zhang Lin Zhu Massimiliano Ferronato Carlo Janna Xulong Gong Jun Yu Pietro Teatini and Shujun Ye

^{1},

^{2},

^{3},

^{4},

^{2},

^{2},

^{3},

^{3},

^{2},

^{1}

^{1}School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China^{2}Dept. of Civil, Environmental and Architectural Engineering, University of Padova, Padova 35121, Italy^{3}Key Laboratory of Earth Fissures Geological Disaster, Ministry of Land and Resources, Geological Survey of Jiangsu Province, Nanjing 210018, China^{4}College of Resource Environment and Tourism, Capital Normal University, Beijing 100048, China

^{1}School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China^{2}Dept. of Civil, Environmental and Architectural Engineering, University of Padova, Padova 35121, Italy^{3}Key Laboratory of Earth Fissures Geological Disaster, Ministry of Land and Resources, Geological Survey of Jiangsu Province, Nanjing 210018, China^{4}College of Resource Environment and Tourism, Capital Normal University, Beijing 100048, China

**Correspondence**: Shujun Ye (sjye@nju.edu.cn)

**Correspondence**: Shujun Ye (sjye@nju.edu.cn)

Earth fissures accompanying anthropogenic land subsidence due to excessive aquifer exploitation create significant geohazards in China. Numerical models represent a unique scientific approach to predict the generation and development of earth fissures. However, the common geomechanical simulators fail to reproduce fissure development because they cannot be effectively applied in discontinuous mechanics. An innovative modelling approach developed recently is applied to develop a software to simulate fissure development. The pressure changes are used as forcing factors in a 3D geomechanical model, which combines Finite Elements and Interface Elements to simulate the deformation of the continuous aquifer system and the generation and sliding/opening of earth fissures. The approach has been applied to simulate the earth fissures at Guangming Village in Wuxi, China with land subsidence of more than 1 m caused by the overexploitation of the second confined aquifer. The modelling results highlight that the earth fissures at Guangming Village have been caused by tension and shear stresses. Based on the developed modelling approach and the application case study, a software platform is developed to provide a fast preliminary evaluation of the risk of fissure occurrence associated to land subsidence. The software allows for the simulation of a simplified 2D conceptual geologic model of earth fissures, which can be used to investigate how the main factors controlling the geomechanical response of the aquifer system, such as pressure changes, geometry of aquifer system, geomechanical properties, and depth of bedrock/fault etc.

Earth fissures accompanying anthropogenic land subsidence due to excessive
aquifer exploitation create significant geo-hazards in China (Ye et al.,
2016, 2018) and worldwide (Conway, 2016; Peng et al., 2016; Teatini et al.,
2018). In China, there are more than 1000 earth fissures identified in
densely urbanized regions, specifically Fenwei Basin, North China Plain, and
Yangtze River delta (Ye et al., 2016). As an example, the direct and
indirect economic loss in Guangming Village of Wuxi city was about USD 20 million in a 1 km^{2} area crossed by two earth fissures (Wang
et al., 2016).

Understanding the generation of earth fissures and modelling their occurrence and propagation are important issues to be addressed. Recently, Ye et al. (2018) proposed a novel modelling approach to simulate earth fissure generation and propagation in 3D complex geological settings. A nested two-scale approach associated with an original non-linear elasto-plastic finite element/interface element simulator allows modeling the mechanics of earth fissures. The approach was successfully applied to the Guangming Village in Wuxi, China, where groundwater pumping since 1980 have caused more than 2 m of land subsidence, and two main earth fissures with maximum sliding and opening of 10 to 40 cm. The model outcomes highlight that the main factors contributing to the earth fissure generation include the shallow bedrock with a sharp ridge, the uneven thickness of the sedimentary deposits at the two sides of the bedrock ridge, piezometric decline and the mechanical parameters of the porous material.

The analysis of the failure processes requires advanced numerical modeling capabilities, and the stable and accurate simulation of earth fissure generation and motion is still an open issue. So it is necessary to develop a software platform to model earth fissures caused by extensive aquifer exploitation based on the approach developed by Ye et al. (2018), which can be easily used by experts working on earth fissures. It can provide a fast preliminary evaluation of the risk of fissure occurrence associated to land subsidence and the main factors that typically control the geomechanical response of the aquifer system, including the geometry of bedrock ridge and aquifer system, pressure changes, mechanical parameters and time scale. The software is developed to simulate simplified 2D conceptual geologic models of earth fissures.

The paper is organized as follows. The numerical models are initially introduced. Then, the developed software and its application are presented. In this study, only the simplified aquifer systems, which represent a conceptual schematization of a real-world hydrogeologic setting at Guangming Village in Wuxi, China, are taken into account. For this kind of conceptual model, earth fissures are generated by tensile fractures above a bedrock ridge (Sheng et al., 2003). A summary section closes the paper.

Based on the pressure changes are known, modelling of earth fissures is carried out with the use of two separate simulations, which are described as: a geomechanical model for the computation of the displacement and stress fields in a continuous domain, and an earth fissure model that addresses the possible generation and propagation of discontinuities within the porous media (Ye et al., 2018).

## 2.1 Geomechanical (GM) model

Based on classical poroelastic theory, the equilibrium equations governing the deformation of a mechanically isotropic medium is:

with *E* the Young modulus, *ν* the Poisson ratio, *u*_{x}, *u*_{y} and
*u*_{z} the displacement components along the coordinate axes *x*, *y*, and *z*,
respectively, ** ε** the volume strain $\mathit{\epsilon}=\frac{\partial {u}_{x}}{\partial x}+\frac{\partial {u}_{y}}{\partial y}+\frac{\partial {u}_{z}}{\partial z}$,

*P*the pressure change and ∇

^{2}the Laplace operator.

The strain ** ε** and stress

**fields can be obtained from the solution**

*σ**u*provided by the GM model using the following relationships:

**D** is the so-called constitutive matrix. The analysis of the stress and
strain fields is important to interpretate the earth fissure development.

## 2.2 Earth fissure (EF) model

A fissure can be considered as an inner boundary embedded in the continuous
body, with a contact condition acting on the opposed surfaces,
Γ^{top} and Γ^{bottom} , that allows for
a relative displacement between corresponding points whenever the stress
state violates the classical Mohr-Coulomb failure criterion:

with *σ*_{n} and *τ* the normal and shear stresses,
respectively, acting on the fissure surfaces, *τ*_{L} the limit shear
strength, *ϕ* and *c* the fissure friction angle and cohesion, respectively.
Any failure criterion prescribes that the perfect continuity of the
displacement across the fissure is preserved if both Eqs. (4) and (5)
are satisfied. Otherwise, the contact surfaces are free to develop relative
movement *u*_{r}:

with *u*_{r} conventionally defined as the movement *u*_{bottom} of
Γ^{bottom} surface with respect to *u*_{top} occurring on
Γ^{top}.

The standard Galerkin Finite Element (FE) method is used to solve the 3D geomechanical models with the presence of active/inactive fissure elements. Fissures are discretized by means of Interface Elements (IEs), with shape functions compatible to those of the surrounding solid elements. The displacements at each time step are calculated using the pore pressure gradient presumed as a known body load. Within each iteration of the non-linear Newton scheme, the resulting 3D displacements are used to compute the stress field and the failure criterion is checked on the IEs: if Eqs. (4) and/or (5) are violated, IE opening and sliding may occur and the fissure develops. More details of the models and numerical methods can refer to the paper by Ye et al. (2018).

## 3.1 Software

The software consisting of the GM and EF models, which is developed for simulating the simplified 2D conceptual cross section models of earth fissures at Guangming Village. The software provides a preliminary evaluation of generation and development of the rupture and the main contributing factors leading to the rupture. The software platform allows for changes in the geologic configuration including the depth of the ridge, and the geometry and thickness of the deposits. It also supports modifying other controlling factors such as pressure, mechanical properties of the lithological units, and time scale. The mathematical models for displacement and stress, and earth fissure development, and their corresponding numerical methods of FE and IE are described in Sect. 2. Fortran (2003) is used to code the models. The software includes modules for Input, Geomechanical properties, Fissure conditions, and Output. All the input files are divided into two folders, INPUT and INPUT_IE. If all the input structure allocation completes and the program runs properly, the outcomes are saved in two folders of OUTPUT and OUTPUT_IE.

## 3.2 Application

To demonstrate the applicability of the software, a simplified conceptual
schematization of hydrogeologic setting in Guangming Village where the
ruptures occur is presented. The domain extends horizontally 2000 m in *X* direction, 50 m in *Y* direction and vertically 500 m in *Z* direction (Fig. 1).
A continuous element mesh is discretized by 288 000 tetrahedral elements and
59 326 nodes. In the software, coarse, medium and fine grids are optional. A
medium grid is chosen in this application.

Displacements are prevented along the direction orthogonal to *x*–*z* plane. So
the software actually simulates the conceptual 2D cross section, although
fissure development is 3D by nature. Roller constraints are prescribed on
the lateral boundaries with a traction-free top plane and zero displacements
on the bottom. The values of mechanical parameters of materials are shown in
Table 1. The aquifer (blue layer) is 350 m thick (Fig. 1). In this layer,
the pore pressure variation is uniform and varies linearly from 0 to −1 Mpa
during a 10 year simulation and recovers inversely to 0 Mpa over the next
10 years. Every year represents a time step, so the whole period contains 20
steps. No pressure change develops in the overlying aquitards (green layer)
and in the sharp bedrock ridge (red area). The contact interface (yellow
line) has 108 triangular elements which are inserted directly above the tip
of the bedrock and reaches the land surface.

The simulation outcomes demonstrate that when the pressure decreases and the presence of the bedrock ridge induces differential subsidence (Fig. 2). The deformation of the compressing sediment above the ridge causes tension (Fig. 3), which is large enough for the generation of the rupture (Fig. 4a). However, the rupture only remains at a shallow depth as its downward propagation is terminated. Figure 4a shows the maximum fissure opening. The vertical displacement of the sediment and the opening both decrease along with the pressure recovery (Figs. 2b, 4b).

This research aims to develop a software platform to numerically simulate earth fissures caused by groundwater withdrawal in a simplified 3D hydrogeological setting of a cross section with a bedrock ridge. The original mathematical models and numerical formulation developed by Franceschini et al. (2016) and Ye et al. (2018) has been used as the mathematical framework for developing the Fortran codes.

The software allows for changes in the geologic configuration of the system including the depth and shape of the bedrock ridge, the geometry and thickness of the deposits, pressure, mechanical properties of the lithological units, and time scale. A case study of a simplified 2D cross section of the hydrogeological setting at Guangming Village, China, demonstrates that the software can be applied to make a preliminary evaluation of the generation and development of a rupture.

The software is still at alpha version. One can contact the corresponding author to get access to the software.

SY, JY, XG and YZ conceived the study and secured funding. SY, YL and PT planned the framework of the software. MF, MF and CJ contributed to the numerical models. YL performed the modeling and develop the software. YL and SY wrote the manuscript. All authors interpreted the results and reviewed the manuscript.

The authors declare that they have no conflict of interest.

This article is part of the special issue “TISOLS: the Tenth International Symposium On Land Subsidence – living with subsidence”. It is a result of the Tenth International Symposium on Land Subsidence, Delft, the Netherlands, 17–21 May 2021.

Funding supported by the Open Project of Key Laboratory of Earth Fissures Geological Disaster of the Ministry of Land and Resources (No. 201703), NSFC (No. 41877180) and UNESCO IGCP641 is appreciated.

This research has been supported by the the Open Project of Key Laboratory of Earth Fissures Geological Disaster of the Ministry of Land and Resources (grant no. 201703), the NSFC (grant no. 41877180), and the UNESCO IGCP (grant no. 641).

Conway, B. D.: Land subsidence and earth fissures in south-central and southern Arizona, USA, Hydrogeol. J., 24, 649–655, 2016.

Franceschini, A., Ferronato, M., Janna, C., and Teatini, P.: A novel Lagrangian approach for a stable numerical simulation of the mechanics of faults, J. Comput. Phys., 314, 503–521, 2016.

Peng, J. B., Sun, X. H., Wang, W., and Sun, G. C.: Characteristics of land subsidence, earth fissures and related disaster chain effects with respect to urban hazards in Xian, China, Environ. Earth Sci., 75, 1190, https://doi.org/10.1007/s12665-016-5928-3, 2016.

Sheng, Z., Helm, D. C., and Li, J.: Mechanisms of earth fissuring caused by groundwater withdrawal, Environ. Eng. Geosci., 9, 313–324, 2003.

Teatini, P., Carreón-Freyre, D., Ochoa-González, G., Ye, S., Galloway, D., and Hernández-Marin, M.: Ground ruptures attributed to groundwater overexploitation damaging Jocotepec city in Jalisco, Mexico: 2016 field excursion of IGCP-641, Episodes, 41, 69–73, 2018.

Wang, G., You, G., Zhu, J., Yu, J., and Li, W. : Earth fissures in Su–Xi–Chang Region, Jiangsu, China, Surv. Geophys., 37, 1095–1116, 2016.

Ye, S., Luo, Y., Wu, J., Yan, X., Wang, H., Jiao, X., and Teatini, P.: Three-dimensional modeling of land subsidence in Shanghai, China, Hydrogeol. J., 24, 695–641, https://doi.org/10.1007/s10040-016-1382-2, 2016.

Ye, S., Franceschini, A., Zhang, Y., Janna, C., Gong, X., Yu, J., and Teatini, P.: A novel approach to model earth fissure caused by extensive aquifer exploitation and its application to the Wuxi case, China, Water Resour. Res., 54, 2249–2269, https://doi.org/10.1002/2017WR021872, 2018.