Abstract
After exploitation of groundwater had been reduced and the groundwater level of the confined aquifer had risen, land subsidence was observed to continue rather than cease for several years according to the layer-wise mark monitoring data in Xi’an. To analyze the phenomena, a numerical model of a coupled one-dimensional multilayer aquifer system is developed to represent land subsidence due to hydraulic head variation in the pumped layer. The numerical simulation results show that the pressure head in other layers does not rise immediately when the hydraulic head in the pumped layer starts to recover after pumping ceases. In addition, after the pumping is stopped, a dividing point can be found in aquitards next to the pumped layer, with the aquitards being divided into two parts: a compressed part and a rebounding part. The dividing points move toward the side and away from the pumped layer with the transferring of pore pressure in the aquitard. The results of the simulation also show that there is a transition period between land subsidence and rebound. In this transition period, land could continue to subside even though the hydraulic head in the pumped layer begins to recover.
1. Introduction
Land subsidence resulting from overexploitation of groundwater has become a major concern all over the world. In China, land subsidence mainly occurs in the 17 provinces in the eastern and central regions, including Xi’an, Shanghai, and the Su-Xi-Chang area [1]. After the exploitation of groundwater had been reduced and the groundwater level of the confined aquifer had risen, land subsidence was observed to continue rather than cease according to the layer-wise mark monitoring data in Xi’an [2, 3]. The same phenomenon also occurred in Shanghai. While groundwater withdrawal in Shanghai had been controlled, the land continued to subside since the 1990s [4]. However, reasons, such as additional stresses caused by building loads and the reduction of recharged groundwater volume, cannot explain the continuous deformation mechanism [5]. It is, therefore, necessary to study the process of land subsidence for analyzing the phenomenon.
The hydraulic head variation in the pumped layer has a significant effect on the process of land subsidence. Therefore, it is necessary to simulate the deformation process of each aquifer in a multilayer aquifer system due to the hydraulic head variation in the pumped layer by applying a numerical model. Several multilayer aquifer system models and numerical analyses have been developed for the analysis of land subsidence. Bear and Corapcioglu [6] presented a variety of analytical solutions to evaluate groundwater flow and solid skeleton deformation due to groundwater withdrawal from aquifer systems. The simulation of multiaquifer systems was carried out using a quasi-three-dimensional model consisting of a three-dimensional hydrologic model coupled with a one-dimensional consolidation model [7–9]. Malama et al. [10] presented a semianalytical solution for the problem of leakage in multiaquifer systems. Considering the effect of the body force and the porosity variation on consolidation, the poroelastic consolidation model of a multilayer aquifer system was developed [11, 12]. The multiscale finite element method (MsFEM) was applied to solve the model of the multilayer aquifer system successfully during the consolidation process [13]. Considering the multiaquifer-aquitard hydrogeological condition, a numerical model was established to analyze relationships among land subsidence, groundwater withdrawal volume, and groundwater level [14]. Lin et al. [15] proposed a distributed model comprised of an analytical quasi-three-dimensional groundwater model and a one-dimensional deformation model to simulate the groundwater flow and deformation for a multilayer aquifer system. Although the study of the process of groundwater flow and deformation in multilayer aquifer systems is extensive, those numerical simulations did not represent the variation state in each layer of the multilayer system during land subsidence. Several experiments of soil consolidation caused by water release have shown that the deformation of the clay layer lags behind the hydraulic head variation. The results from a viscoelastic model developed to simulate one-dimensional consolidation showed that land subsidence obviously lags behind pore pressure disappearance [16]. In addition, the lag time of the clay layer deformation was in direct proportion to the thickness of the layer [17, 18]. The delay index is only dependent on the properties and thickness of the aquitard in response to the sudden change in groundwater level in adjacent confined aquifers according to experimental studies [19].
In this study, a numerical model for a coupled one-dimensional multilayer aquifer system is developed to present the process of land subsidence due to hydraulic head variation in the pumped layer. The variations in pressure head and soil deformation in each layer of the system are simulated to analyze the process of land subsidence along with an analysis of the dividing point that appears in aquitards adjacent to the pumped layer when the hydraulic head in the pumped layer recovers gradually. The movement of the dividing point as a critical part of land subsidence is also explained in this study. There is a transition period between land subsidence and rebound during the hydraulic head recovery in the pumped layer. The transition period could be interpreted as the phenomenon in which land subsidence lags behind the variation in hydraulic head.
2. Methodology
A general mathematical model describing a coupled multilayer aquifer system due to hydraulic head variation in the pumped layer may be formulated by applying Biot’s consolidation theory [20].
The governing differential equations can be derived based on the following basic assumptions: the soil particles are incompressible, the water is slightly compressible, water flow follows saturated Darcy’s law, and the inertial force is neglected in the force equation.
The fluid transport is described by Darcy’s law:where (m/s) is the specific discharge vector, (Pa) is fluid pore pressure, (1/Pas) is the tensor of permeability tensor of the medium, (m2) is the relative mobility coefficient, (kg/m3) is the fluid density, (m/s2) is the gravity vector, and (m) is the position.
The fluid mass balance may be expressed aswhere (1/s) is the volumetric fluid source intensity and is the variation of fluid volume per unit volume of porous material due to diffusive fluid mass transport, as introduced by Biot [21].
The constitutive equation is formulated aswhere (Pa) is Biot modulus, is Biot coefficient, and is the mechanical volumetric strain. The Biot modulus is related to the drained bulk modulus (Pa) and fluid bulk modulus (Pa):where is porosity. If is equal to unity, the grains are considered to be incompressible, and the Biot modulus is equal to.
The relationship between strain rate and velocity gradient iswhere (1/s) is the strain rate. (m/s) and (m/s) are the rate direction of and , respectively.
3. Descriptions of the Multilayer Aquifer System
An ideal seven-layer aquifer system is established to apply the proposed model. The purpose of the numerical model is to simulate the coupling process of pressure head variations and soil deformation in each layer of the multilayer aquifer system subjected to hydraulic head variations in the pumped layer, which is located in confined aquifer II as shown in Figure 1. For such a case, a 300 m thick homogeneous and isotropic multilayer system consisting of confined aquifers I, II, and III, aquitards I, II, and III, and a phreatic aquifer is considered as shown in Figure 1. All the confined aquifers of the model comprise a 40 m thick layer, both aquitard I and aquitard II comprise a 35 m thick layer, the thickness of aquitard III is 25 m, and the phreatic aquifer is 70 m in thickness. The confined aquifers are more permeable and less compressible than the aquitards. In this study, the groundwater of the entire system is assumed to flow only in a vertical direction, and the soil is also assumed to deform vertically.
To describe the relationship between the layer positions and the hydraulic head distinctly, the bottom boundary and land surface of the numerical model are set equal to 0 m and 300 m, respectively. The initial hydraulic head in the saturated zone of the numerical model is set equal to 270 m. It is assumed that the phreatic water table level is constant during all periods so that the vertical downward seepage from the phreatic aquifer recharges each layer directly or indirectly when the hydraulic head in the pumped layer drops. The land surface, also called the top boundary, is free to move vertically. At the bottom boundary, the soil is restrained from vertical movements. Additionally, the no-flow boundary condition is assigned to the bottom side of the model region. The initial displacement and velocity are equal to 0 throughout the model. As a vertical one-dimensional simulation, there is no need to set the lateral boundary. Table 1 lists the hydraulic parameters of each layer implemented, which include the phreatic aquifer, aquitards, and confined aquifers. As mentioned, these parameters are obtained from the references [22] but also are based on average values found in the literature for the Xi’an area undergoing land subsidence.
4. Results and Analyses
4.1. Characteristics and Development of the Land Subsidence Process
A 60-year coupled simulation is distributed in three stages as indicated in Figure 2 based on land surface vertical displacement and hydraulic head in the pumped layer. In addition, the soil is defined as being in compression when the soil vertical displacement is a positive value in this study.
The first stage (Stage I), from the first to the thirty-first year of the entire simulation period, is designed to simulate the fact that the hydraulic head of the pumped layer drops drastically from 270 m to 180 m. The decline of the hydraulic head level induces the land subsidence to begin. In addition, the hydraulic head has been observed to reach the minimum value of 180 m at the end of the first stage. However, the land subsidence did not cease at this stage. The hydraulic head in the pumped layer recovers rapidly in the second stage (Stage II) from the thirty-first to the thirty-eighth year. In contrast to the first stage, the land continues to subside at the beginning of the second stage, and then the land subsidence ceases gradually and reaches a stable state. The vertical displacement reaches the maximum value of 1.16 m at the end of the second stage. The results of the simulation show that the maximum land subsidence lags behind the minimum hydraulic head in the second stage. This indicates that, in the multilayer aquifer system, land subsidence cannot immediately respond to the temporal variation in hydraulic head. The land surface rebounds gradually in the third stage (Stage III) from the thirty-eighth to the sixtieth year. The land surface begins to smoothly rebound at the beginning of the third stage because of aquifer storage and hydraulic head recovery. The hydraulic head recovery occurs smoothly and slowly at the end of the third stage, while the land surface vertical displacement increases rapidly. Notice that the changes in hydraulic head in the pumped layer can only instantaneously influence the deformation of the pumped layer rather than the other aquifers.
4.2. Analysis for Each Aquifer in the Multilayer System
The distribution of pressure head and soil displacement in each layer of the multilayer aquifer system at different stages are shown in Figures 3–8. The head equation, a simplified form of the Bernoulli Principle for incompressible fluids, can be expressed aswhere (m) is the pressure head, (m) is the hydraulic head, and (m) is the elevation head.
(a) Stage I
(b) Stage I
(c) Stage II
(d) Stage II
(e) Stage III
(f) Stage III
(a) Stage I
(b) Stage I
(c) Stage II
(d) Stage II
(e) Stage III
(f) Stage III
(a) Stage I
(b) Stage I
(c) Stage II
(d) Stage II
(e) Stage III
(f) Stage III
(a) Stage I
(b) Stage I
(c) Stage II
(d) Stage II
(e) Stage III
(f) Stage III
(a) Stage I
(b) Stage I
(c) Stage II
(d) Stage II
(e) Stage III
(f) Stage III
(a) Stage I
(b) Stage I
(c) Stage II
(d) Stage II
(e) Stage III
(f) Stage III
The pressure head generally declines and the soil displacement gradually increases in aquitard I during all periods as shown in Figure 3. During Stage I, when the pressure head in the pumped layer (i.e., the confined aquifer II) declines rapidly, the pressure head in aquitard I also declines to recharge the underlying aquifer (see Figure 3(a)). Following the drop in the hydraulic head, the soil displacement in the aquitard I continues to increase in Stage I (see Figure 3(b)). The pressure head in confined aquifer II recovers promptly with groundwater pumping ceasing in Stage II, but none of the pressure head and soil displacement in aquitard I responds to the changes in confined aquifer II (see Figures 3(c) and 3(d)). The pressure head continues to decrease and soil displacement continues to increase with time in aquitard I during Stage III even though the hydraulic head recovery has been maintained for about 30 years in the pumped layer. However, the rates of the pressure head drawdown and soil displacement increase slow (see Figures 3(e) and 3(f)). It is noticed that the variations in the pressure head drawdown and soil displacement initially occur in the lower part of aquitard I. In addition, the soil displacement in the lower part of aquitard I decreases abruptly (see Figures 3(b), 3(d), and 3(f)). The variation in the underlying confined aquifer I could influence first the lower part of aquitard I.
Similarly, the pressure head drops and the soil displacement increases in confined aquifer I at all stages as shown in Figure 4. In other words, the confined aquifer I recharges the underlying aquifer and is in compression all the time. The upper and lower parts of confined aquifer I, which are influenced by the neighboring aquitards, become more compressible. Therefore, the soil displacement in both the upper and lower parts of confined aquifer I is larger than in the middle part (see Figures 4(b), 4(d), and 4(f)).
Aquitard II is the overlying aquifer adjacent to the pumped layer. Large variations in both pressure head and soil displacement initially occur at the bottom of the layer (see Figure 5). The pressure head declines following the hydraulic head drawdown in the pumped layer during Stage I. Meanwhile, the soil displacement increases with respect to the pressure head drop in the aquitard II (see Figures 5(a) and 5(b)). During Stage II, the pressure head decreases gradually with time in the upper part of the layer and increases slightly with time in the lower part of the layer (see Figure 5(c)). In other words, there is a dividing point separating aquifer II into two parts: the upper compressed part and lower rebounding part, with respect to the variation in the pressure head (see Figure 5(d)). It is obvious that the upper compressed part (around 22 meters thick) is thicker than the lower rebounding part (around 18 meters thick) from Figure 5(d). Therefore, aquitard II can be treated as in compression at Stage II. The dividing points of both the pressure head and soil displacement move upward and to the side away from the pumped layer during Stage III (see Figures 5(e) and 5(f)). Because of the upward movement of the dividing point, it is observed from Figure 5(d) that the lower rebounding part (around 30 meters thick) finally becomes thicker than the upper compressed part (around 10 meters thick). Therefore, aquitard II rebounds at Stage III.
As the pumped layer, the confined aquifer II is compressed corresponding to the pressure head drawdown during Stage I (see Figures 6(a) and 6(b)). The upper and lower aquifers recharge confined aquifer II due to pumping being stopped during Stage II. As a result, confined aquifer II rebounds immediately, corresponding to the pressure head rise (see Figure 6(d)). The difference in the pressure head among confined aquifer II and other aquifers becomes smaller with time, and the rate of the pressure head rise and layer rebound becomes smaller.
Aquitard III is the underlying aquifer adjacent to the pumped layer. Large variations in both pressure head and soil displacement initially appear at the top of the layer (see Figure 7). With the pressure head drawdown in the pumped layer, the pressure head declines during Stage I. Meanwhile, aquitard III is compressed (see Figures 7(a) and 7(b)). In contrast to aquitard II, during Stage II, the hydraulic head recovers gradually with time in the upper part of the layer and declines slightly with time in the lower part of the layer (see Figure 7(c)). In other words, there is also a dividing point separating aquifer III into two parts: the upper compressed part and lower rebounding part, with respect to the variation of pressure head (see Figure 7(d)). Aquitard III could be in compression during Stage II because the lower compressed part (around 15 meters thick) is thicker than the upper rebounding part (around 10 meters thick) as shown in Figure 7(d). The dividing points of both the pressure head and soil displacement move downward and to the side away from the pumped layer during Stage III (see Figures 7(e) and 7(f)). The upper rebounding part (around 15 meters thick) finally becomes thicker than the lower compressed part (around 10 meters thick) due to the downward movement of the dividing point so that aquitard III could rebound during Stage III.
As with confined aquifer I, the pressure head declines and the soil displacement increases in confined aquifer III during all stages as shown in Figure 8. Confined aquifer III recharges the overlying aquifer and stays in compression all the time. The upper part of the layer gets more and more compressible influenced by the neighboring aquitard. Therefore, the soil displacement in the upper part of confined aquifer III is larger (see Figures 8(b), 8(d), and 8(f)). Note that the pressure head and soil displacement could not be recorded at elevation of 0 m because the bottom boundary of confined aquifer III is restrained from displacement.
Stage II is a transition period between land subsidence and rebound. As shown in Figure 9, the total compression rate of the entire system decreases rapidly and the total rebound rate increases during Stage II. At the end of Stage II, the total compression rate is equal to the total rebound rate. After that, the land surface starts to rebound. As mentioned above, each layer of the multilayer aquifer system is in compression during Stage I. During Stage II, the performance is different among the aquifers. Confined aquifer II rebounds at once, aquitard I, confined aquifer I, and confined aquifer III are compressed, and both aquitard II and aquitard III are separated into a rebounding part and compressed part by the dividing point. In addition, at the beginning of Stage II, when the total compression deformation is larger than the total rebounding deformation, the land continues to subside even though the hydraulic head in the pumped layer begins to recover. The displacement of the land surface becomes stable (see Figure 2) until the total deformation between compression and rebounding gradually reaches equilibrium at the end of Stage II. The dividing point in aquitard II and aquitard III moves toward the side away from confined aquifer II during Stage III, and the total rebounding deformation is gradually larger than the total compressed deformation. Therefore, the land surface rebounds during Stage III as shown in Figure 2. It appears that the dividing point movement has significant impact on the land subsidence process. The poorly permeable aquitard cannot transfer the pore pressure rapidly relative to the sandy aquifer, which can be observed in Figures 3–8, so that the variation in pressure head could be transferred from one side to the other side in the aquitard for a long period. Therefore, there are different performances between the upper and lower sides of the aquitard. Once one side of the aquitard presents an opposite distribution, the dividing point will appear and move toward the other side with the transference of pore pressure in the aquitard similar to the movement of the dividing points in aquifer II and aquifer III.
5. Conclusions
Pressure head and soil displacement in a seven-layer aquifer system respond to the variation in hydraulic head of the pumped layer during different stages as indicated by a coupled one-dimensional numerical model. The numerical simulation results of the multilayer aquifer system show that the pressure head in other layers will not rise immediately when the pressure head in the pumped layer begins to recover after pumping ceases because the aquitards adjacent to the pumped layer are poorly permeable so that the recovery in pressure head from the pumped layer is transferred to other layers over a long time. The numerical simulation results also show that there will be a dividing point in aquitards adjacent to the pumped layer after pumping ceases. The aquitards adjacent to the pumped layer are divided into a compressed part and rebounding part by the dividing points. The dividing points move toward the side away from the pumped layer over time since the poorly permeable aquitard cannot transfer the pore pressure rapidly relative to the sandy aquifer, and the variation in pressure head can be transferred from one side to the other side of the aquitard for a long period. Therefore, there are different reactions in the upper and lower sides of the aquitard. If one side of the aquitard presents an opposite distribution, the dividing point will appear and move toward the other side with the transference of the pore pressure in the aquitard. The results of the numerical simulation also show that there is a transition period between land subsidence and rebound. In this transition period, land could continue to subside even though the hydraulic head in the pumped layer begins to recover, because, at the beginning of the transition period, the compressed parts are thicker than the rebounding part in the aquitards separated by dividing points adjacent to the pumped layer so that the total compression rate is larger than total rebounding rate in the entire system. At the end of this period, the total compression rate is equal to the total rebounding rate. With the movement of the dividing points, the rebounding parts gradually become thicker than the compressed parts, and the total compression rate is larger than the total rebounding rate. The land surface starts to rebound after the transition period. The conclusions of this study may be helpful in understanding the phenomenon referred to in some studies in which land subsidence lags the hydraulic head in the pumped layer. A fully coupled three-dimensional groundwater flow and land deformation model will be developed to analyze the relationship between the pressure head and soil displacement in future work.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by Xi’an Center of Geological Survey Project “1 : 50,000 Hydrogeological and Engineering geological survey about major district of Guanzhong Basin” (12120113004800).