Modelling ground rupture due to groundwater withdrawal: applications to test cases in China and Mexico

The stress variation induced by aquifer overdraft in sedimentary basins with shallow bedrock may cause rupture in the form of pre-existing fault activation or earth fissure generation. The process is causing major detrimental effects on a many areas in China and Mexico. Ruptures yield discontinuity in both displacement and stress field that classic continuous finite element (FE) models cannot address. Interface finite elements (IE), typically used in contact mechanics, may be of great help and are implemented herein to simulate the fault geomechanical behaviour. Two main approaches, i.e. Penalty and Lagrangian, are developed to enforce the contact condition on the element interface. The incorporation of IE incorporation into a three-dimensional (3-D) FE geomechanical simulator shows that the Lagrangian approach is numerically more robust and stable than the Penalty, thus providing more reliable solutions. Furthermore, the use of a Newton-Raphson scheme to deal with the non-linear elasto-plastic fault behaviour allows for quadratic convergence. The FE – IE model is applied to investigate the likely ground rupture in realistic 3-D geologic settings. The case studies are representative of the City of Wuxi in the Jiangsu Province (China), and of the City of Queretaro, Mexico, where significant land subsidence has been accompanied by the generation of several earth fissures jeopardizing the stability and integrity of the overland structures and infrastructure.


Introduction
The exploitation of subsurface resources, including freshwater aquifers, involves various environmental problems.One of these is land subsidence.Recently, attention has been directed to the activation of pre-existing regional faults, as well as the possible generation of new fractures, caused by groundwater withdrawals.Indeed, the activation of faults can trigger or induce a seismic activity (González et al., 2012).Another problem related to the fault activation is the creation of preferential pathways for fluid leakage.This problem can be also relevant in the case of storage of hydrocarbons or carbon dioxide in the subsoil.Activation of faults and the generation of fissures have a direct effect on the land surface and may cause unacceptable damage to engineered structures and disruption to infrastructures.This is especially problematic in densely populated areas, such as in China (Wang et al., 2007), Mexico (Carreón-Freyre, 2010;Carreón-Freyre and Cerca, 2006), as well as in many semiarid sedimentary basins worldwide (e.g., Azat and Shahram, 2010).A different but somewhat related problem is the hydraulic fracturing or socalled "fracking" technique, which involves the intentional generation of underground fractures in the rock by injecting fluids at high pressure to create preferential pathways that facilitate the extraction of hydrocarbons.A reliable mathematical model to predict the consequences of subsurface fluid withdrawal must take into account the presence of discontinuities, if any, in the porous medium.
At present, the existing numerical-mathematical models to simulate the discontinuities in the porous medium are based on two techniques: (i) use of a continuous model where the fault is characterized by physical properties different from Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.
those of the hosting medium (Rutqvist et al., 2007); (ii) introduction of a discontinuity surface, which behaves as an internal boundary and can be present or absent depending on the stress state.At a numerical level, in the context of the Finite Element Method, the first approach implies a finer discretization of the neighbourhood of the fault, substantially to simulate a layer of different material, while the second requires the introduction of Interface Finite Elements.The first method is simpler, as it involves no real discontinuities, but the difficulties and limitations rest on the proper characterization of a continuous medium that simulates a discontinuity.With Interface Finite Elements (Goodman et al., 1968;Ferronato et al., 2008), there are two approaches to prescribe the necessary mechanical-geometric constraint conditions, that are the non-penetration of solid bodies and the stress state consistent with some failure criterion, e.g., Mohr-Coulomb (Labuz et al., 2012).The first method, referred to as penalty (Bathe, 2006), consists of the introduction of very stiff springs between the faces of the rupture, which impose the condition of non-penetration until they break from overcoming a stress limit.From a mathematical point of view, this method is not exact, because the elastic springs have displacements for any non-zero stress value, while, in nature, in case of a closed fracture the relative displacements are exactly zero.Furthermore, numerically, this method causes a strong ill-conditioning of the stiffness matrix (Ferronato et al., 2012), because of the introduction of the penalty coefficients.However, mainly for the ease of implementation, this approach is widely used.Alternatively the constraint conditions can be imposed by using Lagrange multipliers (Bathe, 2006;Simo et al., 2000), namely in an analytically exact way.The Lagrange multipliers, which physically represent the contact stresses, are additional unknowns, so that the problem is numerically more complex.The increase of the computational cost however, is offset by a more robust convergence in all non-linear steps and a more stable numerical behaviour.
With a Lagrange approach, it is of paramount importance, in case of a rupture sliding, to account for the relative orientation of the shear stress vector with respect to the relative displacement.As the Mohr-Coulomb failure criterion is planar, there is no information on three-dimensionality.A way to compute such an orientation is based on an application of the criterion of "Maximum Plastic Dissipation" (Wriggers, 2006), which uniquely provides the direction of the shear stress as a function of the relative displacement so as to maximize the friction work.By implementing this criterion in the derivation of the variational formulation, we get a numerical scheme that converges quadratically.

Numerical modelling of fault mechanics
A ground rupture is a discontinuity in a porous medium, whose behaviour is governed by the combination of stress and geometry of the close surroundings.To model such behaviour, we choose a failure criterion that prescribes the rupture activation.The shear stress limit τ L is given by (Labuz et al., 2012): where c is the cohesion and σ n the normal stress acting on the fault, taken to be negative in compression.Equation ( 1) defines the modulus but not the direction of the limit shear stress.According to the criterion of "Maximum Plastic Dissipation", the maximization of the friction work W f = τ T L u r , with u r the relative tangential sliding and τ L the limit shear stress, yields: Notice that boldface characters denote vectors.In case of a fault opening, i.e. σ n > 0, the two sides of the rupture become completely free to move with no loading.
To attain equilibrium, the minimization of the functional representing the total potential energy produces the principle of virtual work: where the internal virtual work reads δW i = δε T σ dV , with σ the total stress, and the external virtual work is δW e = δu T bdV + δ δu T tdS, with b and t the volume and surface forces, respectively.The contribution provided by the work of the discontinuity must be added to Eq. ( 3).The total rupture area is denoted by , while the sliding portion, where Eqs. ( 1) and (2) hold, is indicated with ¯ .If opened, the rupture is an unloaded boundary with no contributions to the virtual work.Otherwise, the virtual work δW f reads: As long as the discontinuity is close, both stresses and displacements are independent variables and can have virtual Proc.IAHS, 372, 63-68, 2015 proc-iahs.net/372/63/2015/variations.Conversely, in case of sliding, only the displacement is an independent variable, as the shear stress is fully defined.In this work, we adopt a one-way coupled approach, so the pressure is a priori known and can be treated as an external force, separating it from the term of total stress, according to Terzaghi's principle (Bishop, 1959), that is σ = σ − αpi, where σ are the effective stress, α is the Biot coefficient, p is the pressure and i is the vector form of the Kronecker delta.So we may write the complete conditional formulation as:

Numerical model
The displacements and stresses are approximated by the functions belonging to two different finite Hilbert spaces that, according to the Finite Elements Method, have a dimension equal to the nodal grid number of the entire domain and of the surface of the discontinuity, respectively.In general, we can write: With the assumption of small displacements, deformations are calculated as ε h = Bu h = LNu, where L is the symbolic matrix of differential operators.Denoting by D the tangent stiffness matrix, the variation of the effective stress is dσ h = Ddε h .Moreover, we have to introduce the relative displacements between the two rupture faces, u h r = Su, where the matrix S properly maps the nodes of the global mesh.It is also necessary to express the stress on the discontinuity in a local reference system, using an appropriate rotation matrix R, which binds the global system to the local one.Thus, the discretized form of the Eq. ( 2) becomes: In Eq. ( 7), the normal stress to the rupture surface can be written as σ n = n T λ, where n is the unit vector normal to the discontinuity plane.The functional Eq. ( 5) is now written as: Equation ( 8) must hold for every virtual displacement and stress.So, two non-linear equations are obtained.To solve this system, we use the Newton-Raphson scheme.Computing the Jacobian, we obtain a non-symmetric saddle point matrix.The non-symmetry is caused by the contribution due to the principle of "Maximum Plastic Dissipation" which introduces a component of virtual work where the displacements are independent but the stresses are dependent variables.ments and constant piecewise functions for the Lagrangian multipliers.This choice was performed for compatibility reasons with the stresses of the 3-D FE grid.The discretization relies on the classic tetrahedral FE and the fault is a union of tetrahedra faces.We address two applications: (i) the first concerns earth fissuring in Wuxi, near Shanghai, China, (Wang et al., 2007) caused by groundwater withdrawal from a shallow aquifer, with shallow confined bed; (ii) the second case concerns the City of Queretaro, Mexico (Carreón-Freyre and Cerca, 2006), placed above a graben structure with a volcanoclastic and sedimentary filling (stiff clay, silt and sandy materials) where several earth fissures have been generated by groundwater withdrawal; A sketch of the hydrogeological setting of the two test cases is shown in Fig. 1.For each case, we present the essential geometry and material properties, the history of pumping and the results, with special attention to the behaviour of both subsidence and rupture generation.Preliminary simulations with a simplified stress history are carried out.

Test case 1: Wuxi
In China, earth fissuring related to groundwater pumping has occurred since the 1970s.To cope with the rapid economic development large volumes of groundwater were withdrawal in the Wuxi area, Jiangsu Province, resulting in extensive land subsidence (Shi et al., 2007) and fissure development (Wang et al., 2007).Earth fissures cause serious damages, including cracking of buildings and failure in underground pipelines, with huge economic losses.The geologic setting, which is characterized by an undulating shaped, relatively shallow rocky paleo-basement covered by Quaternary compressible sedimentary deposits from the Yangtze River, strongly enhances the risk of fissure development.
Geological information has been used to develop the static model.A rock ridge, characterized by a depth between few meters and about 100 m from the ground surface is buried below the Quaternary sedimentary sequence (Fig. 2a).The large land subsidence measured in the area between the 1980s and 2000s, has been accompanied by the development of a fissure in correspondence with the tip of the subsurface ridge.The geometry of the main sandy and clayey layers forming the multi-aquifer system has been derived from a number of boreholes drilled down to the bedrock.The target of the study is to simulate the development and propagation of the ground rupture, using a realistic history of groundwater extraction.
The domain extends 2 km × 5 km in the horizontal plane and from the land surface down to 250 m depth in the vertical direction.A refined mesh consisting of 23 303 nodes and 121 942 tetrahedral elements is used to accurately represent the actual lithostratigraphic configuration (Fig. 2b).A num-  ber of 1876 IE are introduced into the 3-D FE grid along the trace of the ridge tip, extending from the land surface to the bedrock top.The IE are characterized by a friction angle φ = 30 • and zero cohesion.The geomechanical properties of the bedrock, aquifers, and aquitards are provided in Table 1.
A uniform pressure decline is prescribed in the sandy layers.The drawdown increases linearly from 0 to 20 m over 10 years.To account for a likely delay of pressure propagation in the clay layers, a half pressure change (i.e., 10 m at the end of the simulation period) is prescribed in the aquitards.No pressure change propagates in the bedrock.
The numerical results confirm the development of an earth fissure along the ridge tip.Due to the quasi-symmetric geological setting, the rupture is characterized by an opening and a negligible sliding.At the end of the simulation period the maximum opening amounts to a couple of decimeters (Fig. 3). Figure 4 provides the corresponding land subsidence distribution.The basement depth of burial significantly affects both the rupture opening and land subsidence.

Test case 2: Queretaro City
Subsidence related groundwater withdrawal and water-level decline has been documented in Queretaro City since the 1970s (Trejo-Moedano and Martinez-Baini, 1991).Local fractures in the fluvio-lacustrine sediments have been reported.The intensity of fracturing has increased and caused numerous problems to urban infrastructure.In the simplified case addressed in this study, an extensive shallow aquifer is covered with a undersaturated and stiff clayey sequence and confined on one side by a hard rock mass.The rock mass is separated from the aquifer by a regional pre-existing vertical fault.The clayey sediments have a thickness of 20 m, while the aquifer extends for 30 m in depth.The mesh covers a 30 km-side square area, one third of which is occupied by rock.In the z direction, the mesh extends from the free surface down to 50 m.Pumping takes place from the aquifer through three wells, as shown in Fig. 5, yielding a maximum pressure drawdown of about 20 m.In 10 years, a linearly increasing flow rate is implemented, up to the value of 300 L s −1 per well.The geomechanical parameters are provided in Table 1.Permeability of clayey sediments and rock is K = 10 −10 m s −1 , while permeability of aquifer is K = 10 −4 m s −1 .The fault has zero cohesion and friction angle φ = 30 • .The FE mesh consists of and 297 971 tetrahedra with 57 050 nodes.The fault is discretized by 8074 IE.
The model results provided by the IE and FE are shown in Figs. 6 and 7, respectively.In this geological setting the fault mainly slides.It can be noted that the fault opens negligibly (less than 3 cm), while the sliding in the fault plane approaches 40 cm, a value of the same order of magnitude as the maximum recorded subsidence, as can be seen in Fig. 7.This shows how important the presence of the fault is for the exact prediction of the land subsidence.The greatest calculated displacements along the fault are in the upper central part, i.e., in areas where the influence of the wells is largest and compaction of the underlying layers is the greatest.The displacements are quite regular and representative of the real order of magnitude.

Conclusions
The main subject of this investigation is how to simulate a geological discontinuity in a continuous medium.We used the Interface Finite Element method with a Lagrangian approc-iahs.net/372/63/2015/Proc.IAHS, 372, 63-68, 2015 proach.At variance with existing algorithms, we introduce the principle of "Maximum Dissipation Plastic" in the formulation of the functional to be minimized.The solution is generally very smooth, as can be seen from the results discussed above.
The algorithm was implemented and preliminarly tested in two cases: (i) an aquifer with a shallow undulating bedrock representative of the geological setting in Wuxi, China.and (ii) a laterally rock confined aquifer covered with stiff clayey sediments, representative of the typical geological setting in Queretaro City In both cases the algorithm is robust and provides numerically stable and physically plausible solutions.Further improvements will be carried out for both the case studies, with the implementation of a more realistic evolution of the pressure head the calibration of FE-IE models on available measurements, and the development of future scenarios accounting for the expected management of the subsurface water resources.

Figure 1 .
Figure 1.Sketches of the typical geological setting ground rupture in Queretaro City (a) and Wuxi, China (b).The simulated configurations are highlighted by the red boxes.

Figure 2 .
Figure 2. Test case 1: perspective views of the geological setting (a) and the 3-D FE-IE grid (b).The materials are highlighted in various colours: rock in red, aquifers in blue, and aquitards in green.

Figure 3 .
Figure 3. Test case 1: opening and sliding of the ground rupture after 10 years of water withdrawal.

Figure 4 .
Figure 4. Test case 1: land subsidence after 10 years of water withdrawal.

Figure 5 .
Figure 5. Test case 2: 3-D mesh with the materials highlighted in various colours.Rock is red, aquifer is blue, and clay is green.The fault divides rock from the sedimentary clay-aquifer system.

Figure 6 .
Figure 6.Test case 2: opening and sliding on the fault after 10 years of water extraction.

Figure 7 .
Figure 7. Test case 2: land subsidence after 10 years of water extraction.

Table 1 .
Geomechanical properties used for the two test cases.E: Young modulus; ν: Poisson's ratio; e: void ratio.