Land subsidence caused by a single water extraction well and rapid water infiltration

Abstract. Nowadays several parts of the world suffer from land subsidence. This setting of the earth surface occurs due to different factors such as earth quakes, mining activities, and gas, oil and water withdrawal. This research presents a numerical study of the influence of land subsidence caused by a single water extraction well and rapid water infiltration into structural soil discontinuities. The numerical simulation of the infiltration was based on a two-phase flow-model for porous media, and for the deformation a Mohr–Coulomb model was used. A two-layered system with a fault zone is presented. First a single water extraction well is simulated producing a cone-shaped (conical) water level depletion, which can cause land subsidence. Land Subsidence can be further increased if a hydrological barrier as a result of a discontinuity, exists. After water extraction a water column is applied on the top boundary for one hours in order to represent a strong storm which produces rapid water infiltration through the discontinuity as well as soil deformation. Both events are analysed and compared in order to characterize deformation of both elements and to get a better understanding of the land subsidence and new fracture formations.


Introduction
Land subsidence processes can occur due to extraction of minerals in mining galleries, tunnels construction, fluids extraction (water, oil or gas) from natural reserves, decrease of groundwater level during prolonged extraction, natural land dissolution, compaction of soil materials or tectonic activity.In modern times, cities have grown rapidly and likewise have the agricultural and industrial activities.In order to satisfy the growing demands for the vital liquid, for which the superficial water bodies have not been sufficient, groundwater has been exploited to such an extent that its extraction has outpaced groundwater recharge.Especially in arid and semi-arid zones this has produced a deficit leading to a rapid decline of the groundwater table.In heterogeneous strata the fast decline of the groundwater table leads to a differential adjustment and in some cases to a formation of cracks and fractures on the surface of the ground.The valley of Queretaro city, Mexico (Fig. 1, right), is an example of active land subsidence.Since the 1970s, with the rapid growth of the city, the water demand increased rapidly and as a result the urban area is affected by differential compaction and formation of a reticular system of faults and fractures.Many of them appear on the surface and have caused economical damages in the last 40 years.
Groundwater overexploitation causing water table decline in unconfined aquifers leads to deformation of the porous media.This is mainly controlled by the increase in effective stress in the soil mass, reducing pore volume and, hence, the aquifer will compresse causing subsidence (Rojas et al., 2002;Martinez et al., 2013) and may cause the generation of earth fissures at land surface.Fault and fractures can influence groundwater flow in aquifer and aquitard layers (Carreón-Freyre et al., 2005).Structural discontinuities can perform as hydrological barriers.Such a barrier can result in different water levels on either side of the fault and as a result differential land subsidence occurs across the fault on the surface of the earth.Furthermore, it has been shown that fast water infiltration through a pre-existing fracture zone into a soil system could perform as a mechanism for fracturing and triggering of land subsidence (Martinez et al., 2013).Groundwater depletion occurs at scales ranging from a single well to regional aquifer systems.The extents of the resulting effects depend on several factors including pumpage and natural dis-Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.charge rates, physical properties of the aquifer, and natural and human-induced recharge rates.
Figure 1 left shows the urban area of Queretaro city (light green), faults and fractures (red line) and some of the active wells in the city (circles).In many of the cases the wells at both sides of the faults show a considerable difference of piezometric levels.An example of different water levels across a foult are given by, the static levels of the wells number 612-F at the left side of the fault Tlacote and the well number 1313-A of the year 2011 are 153.9 and 140 m respectively, which is a difference almost of 14 mm.The same behaviour exists also in wells number 1322 and 1638 located at the right and left sides of the fault Epigmenio Gonzales, respectively, where the water depths in 2011 were 116 and 129.5 m, respectively (Table 1).
The hydrogeological system is complex because of the heterogeneity of the soil layers, the basement, the fracture system and the well system.Beside the complexity of the system, another problem is the lack of data.In the literature there are several examples of analysis of subsidence produced by one single well, for example Ziaie and Rahnama (2007) presented in a study the relationship between classical soil parameters and subsidence driven for a pumping single well.Also Martinez et al. (2013) demostrated that fast water infiltration could be a mechanism for fracture generation and for triggering land subsidence.In this research both elements are analysed.
The aim of this research is to simulate the dynamic groundwater table of a homogeneous aquifer with a fault zone which has been caused by a single water extraction well and fast water infiltration after a strong rainfall.Although only idealized hydrogeological systems are modelled, the re- sults help us to understand the land subsidence intensity resulting from both events, extraction and rapid infiltration and the interrelation of them.

Two-phase flow in porous media
If two fluid phases are not or only slightly miscible into each other, a two-phase flow model concept for porous medium can be applied (Hinkelmann, 2005).The continuity equation must be fulfilled for each phase α, one for the liquid phase w (water) and one for the gas phase n (air): Proc.IAHS, 372, 33-38, 2015 proc-iahs.net/372/33/2015/Here, ϕ α is the void space filled with phase α, S α the saturation, ρ α the density, v α the filter or Darcy velocity and q α a sink or source of the phase.
To describe the Darcy velocity of each phase, the generalized form of the Darcy law can be used: where k rα is the relative permeability, K the tensor of intrinsic permeability, µ α the dynamic viscosity, p α the pressure and g the gravity.k rα µ α = λ α is referred to as the mobility of phase α.
The relative permeability for liquid phase w and gas n phase are computed with the Brooks-Corey (Brooks and Corey, 1964) relationship: The parameter λ characterizes the grain-size distribution.A small value describes a single grain-size material, while a . (5) The pore volume is completely filled with the wetting and the non-wetting phase saturation, S w , S n : At the interface between both, the difference between the phase pressure of gas p n and liquid phase p w is called capillary pressure.
The authors would like to mention that a Richards model concept also would have been suitable for the simulations.Further information about the model is found e.g. in Hinkelmann (2005

Linear-elastic perfectly plastic model with Mohr-Coulomb failure criterion: deformation
Plasticity refers to irreversible strains.As soil behaviour is highly non-linear and irreversible, a linear-elastic perfectly plastic Mohr-Coulomb model has been applied (Waterman et al., 2004).This model involves five input parameters, i.e.Young's modulus E, Poisson's ratio υ, friction angle ϕ, dilatancy angle ψ and cohesion C. For simplification, in this investigation the dilatancy angle is taken equal to zero.Material models for soil and rock are generally expressed as a relationship between infinitesimal increments of effective stress ("effective stress rates") and infinitesimal increments of strains ("strain rates").This relationship is expressed in Hooke's law: where M is a material stiffness matrix.This model decomposed strains and strain rates into elastic and plastic parts: If Eq. ( 8) is inserted into the Hook's law Eq.( 8) leads to: Ḋe is the elastic material stiffness matrix σ effective stress rates.
A plastic potential function g is introduced.The plastic strain rates are written as: in which λ is the plastic multiplier.For elastic behaviour λ is zero and for plastic behaviour λ is positive.These model concepts are implemented in the tool PLAXIS.Further information about the model is found in Waterman et al. (2004).

Idealized system and parameter
A rectangular system of 40 m length and 30 m depth was simulated.A 2-D steady-state model of a two-layered system with a clay layer on the top, a sand layer underneath and a 1 m fault zone was set up.A fracture zone was set up with a 45 • inclination with respect to the horizontal (Fig. 2, left).In this study the fault zone was idealized as a porous medium damage band with higher permeability (blue in Fig. 2, middle) compared to the sand and the clay layers (green and red, respectively, in Fig. 2, middle).An unstructured grid was generated with 9798 cells with a higher resolution around the fault zone as shown in Fig. 2 middle.The boundary condition (BC) at the top of the system is a Dirichlet BC with a fixed atmospheric pressure.At the bottom, left and right of the domain a fixed Neumann no flow BC was chosen.For all simulations, as an initial condition, full water saturation was considered.
The flow parameters used in the numerical simulation for the sand, clay and fracture zone are described in Table 2.

Reference case
Firstly a system without fracture was modelled as a reference to be compare later with the results of the other setups.Figure 3 shows the normal water depletion cone that appears around a well extraction.

Water withdrawal of a single well
With difference to the reference case, in the first step of simulation, the water extraction of a single well near a high permeable fault zone, a barrier effect is produced, i.e. the water level at both sides of the fault zone are different.The results of water saturation of the system with a 45 • high permeability fault zone and a single well extraction (Fig. 4) show a difference of approximately 3.5 m on both sides of the fault zone.This barrier produces an irregular soil deformation, especially on the top of the surface (Fig. 5, left).The results of the simulation show a maximum total vertical displacement of approximately 1.8 × 10 −3 m (Fig. 5, right) and approximately a 0.6 × 10 −3 m maximum total horizontal displacement (Fig. 6).

Fast water infiltration
To simulate an accumulation of water caused by a runoff after a strong rainfall a 1 cm water column on the complete surface was applied.The BC at the top of the system is a Dirichtlet BC with a fixed water pressure p n = p atm + ρgh with h = 1 cm simulating an accumulation of a 1 cm water column (Fig. 2, right).
As shown in Fig. 7, which highlights the deformed mesh after fast water infiltration, the maximal displacements take place in the upper part of the fault zone.
The results show that after 30 min of fast water infiltration through the high permeable fault zone a maximum horizontal displacement of approximately 1.8 × 10 −3 m and a vertical displacement of 6.8 × 10 −3 m take place (Fig. 8).

Conclusions
The present study is a conceptual model of land subsidence produced water extraction from a single well and fast water infiltration through a fault zone.
A two-layered system consisting of a sand layer and a small clay layer on the top with an inclined high permeable fault zone was investigated.
The differences of the water levels on both sides of the fault zone produced by the single well water extraction lead to different settlements on both sides.However, there is a greater impact on the vertical than on the horizontal displacement.The fast water infiltration impacts the displacements in both directions, vertical and horizontal, and the horizontal displacement can cause plastic deformation and also could produce new fracture formations.
If the fault zone is narrower, for example 40 cm, the infiltration through the fault zone is not significant and the water spreads out more into the partially saturated zone because of capillarity forces and gravity (Martinez et al., 2013).This phenomenon increased the deformation of the soil layers, specially the horizontal one.The results also show that a fault zone with high permeability could act as hydrological barrier.This phenomenon will be studied in more detail in future work.

Figure 1 .
Figure 1.Queretaro, urban area, well distribution and faults and fractures location after Pacheco (2007); State of Querétaro within Mexico.

Figure 2 .
Figure 2. Idealised system (left).Unstructured grid of a 40 m×30 m system with higher resolution around the fault zone (middle); Complete 1 cm water column on the top boundary (right).

Figure 3 .
Figure 3. Water saturation after water extraction of a single system without fracture.

Figure 4 .
Figure 4. Water saturation in the case of a single well extraction with a fault zone.

Figure 6 .
Figure 6.Maximum total horizontal displacement after water extraction.

Figure 8 .
Figure 8. Maximum total vertical displacement after fast water infiltration (left).Maximal total horizontal displacement after fast water infiltration (right).

Table 1 .
Water depths of some wells in Queretaro in 2011.

Table 2 .
Soil parameters for two-phase flow modeling.