Next Article in Journal
Hydrogeochemical Characteristics and Health Risk Assessment of Groundwater in Grassland Watersheds of Cold and Arid Regions in Xilinhot, China
Next Article in Special Issue
Morbidity and Water Quality: A Review with a Case Study in Tonosí, Panama
Previous Article in Journal
Effects of Formulated Pellet Feed or Live Fish Food on the Intestinal and Aquaculture Water Microbial Communities in Goldfish, Carassius auratus
Previous Article in Special Issue
Groundwater Sustainability Planning in California: Recommendations for Strengthening the Kern Groundwater Sustainability Plan
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Prediction of Ground Subsidence Induced by Groundwater Mining Using Three-Dimensional Variable-Parameter Fully Coupled Simulation

1
Nanjing Geological Survey Center, China Geological Survey, Nanjing 210016, China
2
Geological Survey of Jiangsu Province, Nanjing 210018, China
3
Key Laboratory of Earth Fissures Geological Disaster, Ministry of Natural Resources, Nanjing 210000, China
4
School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China
5
School of Earth Science and Engineering, Sun Yat-sen University, Zhuhai 519082, China
*
Author to whom correspondence should be addressed.
Water 2024, 16(17), 2487; https://doi.org/10.3390/w16172487
Submission received: 2 August 2024 / Revised: 29 August 2024 / Accepted: 29 August 2024 / Published: 1 September 2024
(This article belongs to the Special Issue Studies on Water Resource and Environmental Policies)

Abstract

:
In order to predict the ground settlement in a scientific, intuitive, and simple way, based on the theory of Bio-consolidation, a three-dimensional fluid-solid coupled numerical calculation programme FGS-3D for ground settlement was compiled by using the Fortran 95 language, and a front-end operation platform was developed by using Microsoft VisualBasic, so that a three-dimensional variable-parameter fully coupled viscoelastic-plastic model of ground settlement was constructed using the city of Yancheng as an example, and the development of ground settlement and horizontal displacement changes from 2021 to 2030 were predicted. The results show that the three-dimensional fully coupled finite-element numerical model of building load, groundwater seepage, and soil deformation established by the above computer development program can directly create a hydrogeological conceptual model of groundwater mining and predict ground settlement, so as to achieve the visualisation of the three-dimensional seepage of groundwater and the fully coupled simulation of ground subsidence in the whole process of groundwater mining. Under the joint action of construction load and groundwater mining, the water level of the aquifer in Yancheng City rises by 1.26 m on average in the main groundwater mining area of the group III pressurised aquifer, forming two smaller landing funnels, and the lowest water level of the two landing funnels is −15 m.

1. Introduction

Ground subsidence hazards are concerns of all corners of society, including the economy, law, groundwater resources, flood control planning, and urban construction, and cause significant changes in urban geomorphological patterns and failure of the original level elevation. In a way, controlling ground subsidence is not only a geological disaster problem, but also a regional environmental problem [1].
Prediction methods for ground settlement caused by groundwater extraction can be divided into empirical, semi-theoretical, and theoretical methods [2,3,4,5,6,7,8,9,10,11]. Ground settlement deformation caused by the release of water from a weakly permeable layer is difficult to portray without establishing water flow equations and modelling water level changes. At the same time, if the coupling of groundwater seepage and ground settlement is not realized mechanistically, only the groundwater seepage model is coupled with the ground settlement model through a certain physical quantity, which is only a partially coupled model, especially if the three-dimensional nonlinear deformation of the soil layer is not taken into account [12], and the soil deformation models adopted are therefore all one-dimensional linear elastic models, which are not able to simulate the prediction of horizontal deformations of the soil layer [13,14].
Currently, the mathematical models used to simulate the prediction of groundwater mining and ground subsidence, both domestically and abroad, are mainly based on Darcy’s law and Terzaghi’s one-dimensional consolidation theory, which are based on the coupled model of groundwater mining and ground subsidence (Romolo Di Francesco, 2013 [15]; Rogers, 2016 [16]), in which the coupled model of groundwater mining and ground subsidence is the most important one. Xu et al. proposed a numerical model which considered 3D seepage and 1D consolidation [17], but the numerical model could not predict the horizontal displacement. This method failed to achieve full coupling. Zhang et al. predicted the land subsidence using Biot’s theory [18], but the constitutive model only considers elastic deformation. Wang et al. constructed a 3D model using TOUGH and FLAC3D equipped with a thermo-elastoplastic constitutive model [19]. The model was used to reproduce groundwater flow, heat transfer, and mechanical responses in porous media. However, the numerical model did not consider the variation of parameters.
In this study, we presented a three-dimensional coupled numerical analysis programme for groundwater mining and ground subsidence, FGS-3D, in Fortran 95. The parameters of the soil are not treated as constant in the programme. Yancheng City, a region located in China, is selected as an example. The groundwater level and the soil deformation are calculated by the programme [20,21,22].

2. Materials and Methods

In this study, a three-dimensional fully coupled finite-element program FGS-3D for groundwater mining and ground subsidence was developed in Fortran 95, a software system that has added new procedures for the overall manipulation of arrays, internal functions of arrays, and dynamic arrays, with the newest feature being the ability to perform modular operations. A modular operation is a programme unit but separate from the main programme, and contains a set of subroutines and functions constituting a “function library” in order to be readable and portable by the operator, just like the difference between subroutines and functions and the main programme.

2.1. Biot’s Consolidation Theory

The three-dimensional water-soil fully coupled mathematical model, as shown in Equation (1), is established based on the more rigorous theory of Biot’s consolidation [23,24,25,26,27,28].
G 2 w x G 1 2 ν x w x x + w y y + w z z + u x = 0 G 2 w y G 1 2 ν y w x x + w y y + w z z + u y = 0 G 2 w z G 1 2 ν z w x x + w y y + w z z + u z = γ t w x x + w y y + w z z + 1 γ w x k x x u x + y k y y u y + z k z z u z + γ w + W = 0
where G is the shear modulus; ν is the Poisson ratio; w x , w y , and w z are the displacement components in the x, y, and z directions; u is the pore water pressure; k x , k y , and k z are the hydraulic conductivities in the x, y and z directions; γ is the specific weight of soil; γ w is the specific weight of water, and 2 is the Laplace operator.
Equation (1) can be abbreviated as:
B T D B w + B T M u = f C T t w 1 r w C T K C u = Q
where w = w x , w y , w z T : x , y , z are the displacements of the computational point in the corresponding direction; D is the matrix of stress-strain relationships (the constitutive equations); f = f x , f y , f z T : x , y , z is the volume force at that point in the direction of the body force; u is the hyperstatic pore water pressure; M = 1 , 1 , 1 , 0 , 0 , 0 T ; the permeability coefficients are k x , k y and k z in three separate directions; Q is the equivalent flow rate at the point of calculation; C = N i x N i y N i z T i = 1 , 2 , , 8 ; γ w is the heaviness of water; K = k x 0 0 0 k y 0 0 0 k z ; and B i = N i z 0 N i y 0 0 N i x 0 N i z N i x 0 N i y 0 N i x N i y 0 N i z 0 0 T i = 1 , 2 , , 8 .
The soil constitutive relationship is a mathematical expression describing the stress-strain relationship of the soil. Rheologically characterised soils are viscoelastic-plastic bodies with elastic, plastic, and viscous properties. The total strain increment d ε is the elastic d ε e , viscoelastic d ε v e , and viscoplastic strain increments d ε v p . Then, the strain increment when considering the rheological properties of the soil can be expressed as:
d ε v = d ε e + d ε v e + d ε v p
1. Dynamic modelling of parameters
(1)
Non-linear relationship between the permeability coefficient and effective porosity
Under the influence of construction loading and groundwater mining, pore water stress dissipation causes soil skeleton stress redistribution, and the stress state in the soil changes, which is macroscopically manifested in the consolidation and deformation of the soil, and in changes in parameters such as effective porosity n and permeability. According to the Kozeny–Carman equation, combined with the theory of Bio-consolidation and the actual situation in the study area, ignoring the changes in the surface area of the geotechnical particles and the small changes in the volume of the thermally expanding particles, we can obtain the equation of the relationship between the permeability, the body strain, and the initial effective porosity as follows:
k = k 0 1 + ε v 1 + ε v n 0 3
where n 0 is the initial porosity, k 0 is initial penetration rate. and ε v is body strain.
(2)
Nonlinearity of the deformation modulus and Poisson’s ratio
According to Duncan–Chang model, the expressions for tangential deformation modulus and tangential Poisson’s ratio are as follows:
E t = 1 R f 1 sin φ σ 1 σ 3 2 c cos φ + 2 σ 3 sin φ 2 k p a σ 3 p a m
ν t = β F lg σ 3 p a 1 D σ 1 σ 3 α p a σ 3 p a m 1 R f 1 sin φ σ 1 σ 3 2 c cos φ + 2 σ 3 sin φ 2
2. Fixed Solution Condition
(1)
Initial conditions
① The initial condition of the pore water pressure is the flow field of each aquifer on 1 January 2016, the value of which is given by the actual measurement, and the initial flow field of the weak aquifer between each aquifer is given by the interpolation of the upper and lower aquifers:
u ( x , y , z , t ) t = 0 = u 0 ( x , y , z ) = γ w h 0 ( x , y , z , t 0 ) ( x , y , z ) Ω
where u 0 ( x , y , z ) is the known initial pore water pressure in the study area, h 0 ( x , y , z , t 0 ) is the initial head value, γ w is the heavy for water, and Ω is the entire study area.
② Geostress initial conditions
The self-gravity stress of the soil was used to estimate the initial stress of the soil:
σ z = γ z σ x = K 0 γ z
where σ x and σ z are the initial horizontal and vertical stresses in the soil body, z is the depth of the calculation point for each layer, K0 is the static lateral pressure coefficient K 0 = 1 sin ϕ 0.95 sin ϕ sand clay , and φ is the effective angle of internal friction.
③ Displacement initial conditions
The initial displacement values for each aquifer on 1 January 2016 were set to 0 as the initial condition for displacement in the model:
w ( x , y , z , t ) t = 0 = 0           x , y , z Ω
(2)
Boundary conditions
① Flow boundary conditions Γ 2 are as follows:
K H n   Γ 2 = q L ¯
where q L ¯ is the known flow rate per unit area Γ 2 on the boundary.
② Free surface boundary conditions Γ 3 are as follows:
u = Z ;         q = μ u t cos θ
where μ is the degree of soil watering, θ is the angle of intersection of the direction of the normal to the outside of the free surface and the vertical line, q is the flow rate per unit area through the free surface boundary Γ 3 , and Z is the elevation at which the free surface is located.
③ Displacement boundary conditions Γ 4 are as follows:
w x Γ 4 = w x ¯ w y Γ 4 = w y ¯ w z Γ 4 = w z ¯
where w x ¯ , w y ¯ , w z ¯ are the known displacements in the three directions on the displacement boundary Γ 4 , and the displacements in the three directions on the displacement boundary are taken to be 0 for this simulation.

2.2. Research on the Development of the Software Programme

The three-dimensional fully coupled numerical analysis program of the building (structural) loading, groundwater mining and ground settlement was compiled using Fortran 95 [29,30,31,32], of which the structure block diagram of the free surface iteration algorithm is shown in Figure 1, that of the viscoelastic-plastic stress analysis iterative algorithm is shown in Figure 2, and that of the fully coupled analysis procedure is shown in Figure 3.
In order to make the above groundwater unsteady seepage three-dimensional finite-element computer programme and the groundwater seepage and soil deformation fully coupled three-dimensional finite-element computer programme have better application value, stronger operability, and better expression of the calculation results, we developed the groundwater system model (GSM) using the Microsoft Windows 11 operating system environment and Microsoft VisualBasic 6.0 language The software consists of seven components: file, modelling, data module, calculation module, post-processing module, restore, and help. The functions of each constituent module are shown in Figure 4.
The drop-down menu includes options for both AutoCAD models and ANSYS models for pre- and post-processing graphic files. The modelling module generates 3D stereo modelling and raw calculation files, mainly including identification of the model, rectangular area units and node information, initial and boundary conditions, and well information generation. The data module loads and modifies model-specific parameters, including parameters, initial conditions, boundary conditions, well files, solutions, and time steps. The calculation module selects different calculation methods to calculate and solve the finite element model in different ways. The drop-down menu includes two items: a three-dimensional seepage calculation and three-dimensional Beal consolidation seepage and stress field fully coupled calculation, compiled by the Visual Fortran 95 language. In this study, the fully coupled calculation method of groundwater extraction and ground subsidence is used. The post-processing module is mainly for providing output and graphing of calculation results, including generating dynamic curves, outputting data, and drawing contour plots.
The hydrogeological conceptual model of groundwater mining can be created directly on the computer, and identification of the model and forecasting of the model can be carried out so as to achieve the visualisation of the fully coupled simulation and calculation of the three-dimensional seepage of groundwater and ground settlement in the whole process of groundwater mining.

3. Model Setup

Based on the visual computer model (GSM) developed in this study, in order to verify its correctness and practicability, it was practically applied to the Yancheng City Ground Subsidence Project.
Geological profile: Yancheng City is located in the eastern coastal open area of Jiangsu Province, in the middle of the North Jiangsu Plain, east of the Yellow Sea, with rich land, sea, and mudflat resources, and is the prefecture-level city with the largest land area and the longest coastline in Jiangsu Province. The study area is mostly a vast plain area with no bedrock outcrops, dominated by loose rock-like pore groundwater, which is very deeply buried, mostly greater than 1000 m. Details of the division of the aquifer groups are shown in Table 1.
Discretisation of the model: According to the structural characteristics of the Quaternary stratigraphy in Yancheng City, the study area was dissected using a three-dimensional model. The entire study area can be divided into 4355 rectangular grid cells in the plane, and vertically into a total of seven separate layers from top to bottom for the calculations, including submersible aquifers, group I pressurised aquifers, group II pressurised aquifers, group III pressurised aquifers, and weak aquifers of clayey soils between the aquifers. The vertical depth ranges from 230 to 350 m and covers an area of 2130 km2. The total number of units is 30,485 and the total number of nodes is 36,984 (Figure 5).
Identification of model: Using the finite-element computer model developed here, observations from 1 January 2016 to 31 December 2020 were selected for inverse analysis of the parameters of the model. Every 3 months is a stress cycle, and the five years were divided into a total of 20 stress periods, with each stress period being one time step. After identification of the model, the model was dissected into 54 parameter partitions. Taking the group III pressurised aquifer as an example, the parameter zoning is shown in Figure 6 and the parameter values are summarised in Table 2. The calculated flow field in the confined group III aquifer as of 31 December 2016 is shown in Figure 7, and the calculated flow field as of 31 December 2018 is shown in Figure 8. The calculated amount of ground settlement over the two years from 1 January 2016 to 31 December 2018 is shown in Figure 9. The fit for the 2016–2018 levelled observations is shown in Figure 10.
According to our previous study, the deformation modulus (E), cohesion (c), and specific weight of soil ( γ ) have significant impacts on land subsidence, and the vertical hydraulic conductivity ( K z z ) has a minimal impact on land subsidence [33].
The water balance elements per year for the group III confined aquifer from the model identification calculations are detailed in Table 3. From the table, it can be seen that the groundwater inflow and outflow of each aquifer group in each regime period are basically equal, which indicates that the model has good convergence and stability, and it can be used to predict the characteristics of the stress-strain dynamics of the groundwater system in the future.

4. Results and Discussion

Prediction of the model: The ground subsidence model was used to simulate and predict the development trend of the groundwater flow field and ground subsidence in various aquifers in Yancheng City from 1 January 2021 to 31 December 2030 under the joint action of building loads and current groundwater mining.
The average rise in the water level of the group III pressurised aquifer from 1 January 2021 to 31 December 2030 is 1.26 m when considering the effects of building loads and current groundwater extraction. In the main groundwater extraction area of the group III pressurised aquifer, two smaller landfall funnels are formed: first is the area bounded by the towns of Qinnan and Xuefu, and second is in the vicinity of Lungang Township. The minimum water level in these two landing funnels is −15 m. However, the group III pressurised aquifer still shows a regional large funnel centred on Longgang Town–Guo Mang Town–Qinnan Town–Xuefu Town. The predicted cumulative maximum ground rebound in Yancheng City from 1 January 2021 to 31 December 2030 is 81.5 mm, and the maximum ground rebound rate is 8.15 mm/a, which is mainly located in Dagang Town and Bicang Town. The predicted cumulative maximum ground settlement is 3 mm, and the maximum ground settlement rate is 0.3 mm/a, mainly located in the vicinity of Chengbei District and Longgang Town. Due to the small amount of extraction, the groundwater level in the study area has continued to rebound, and rebound has occurred in the study area as a whole, but in the centre of the city, where the building loads are high, the amount of ground rebound has been small, and even slight subsidence has occurred.
Details of the predicted flow field map, compressibility contour map, and horizontal displacement map of the group III pressurised aquifer are shown in Figure 11, Figure 12 and Figure 13; the 10-year cumulative ground subsidence contour map is shown in Figure 14.
Discussion: When simulating the change in groundwater level from 1 January 2021 to 31 December 2030, an analysis of Figure 7 and Figure 8 showed that the water level of the group III pressurized aquifer rises by an average of 1.26 m. In the main groundwater extraction area of the group III pressurized aquifer, two smaller landfall funnels were formed: the first one is the area bounded by Qinnan Town-Xuefu Town, and the second one is in the vicinity of Longgang Town, where the lowest water level is −15 m. In addition, a regional large funnel centred on Qinnan Town-Guo Mang Town-Xuefu Town is still shown. In addition, the group III pressurised aquifer still shows a regional large funnel centred on Longgang Town–Guo Meng Town–Qinnan Town–Xuefu Town.
From the contour map of the compression volume from 1 January 2021 to 31 December 2030, it can be seen that the maximum compression volume of the Yancheng Submerged Aquifer is 14 mm, located in the west area of the city; the maximum horizontal displacement is 3 mm, mainly located in the west area of the city and the south area of the city. The group III pressurised aquifer undergoes a relatively large rebound, with a maximum rebound of 16.5 mm, mainly located in Xinxing Town, and a minimum rebound of 5.5 mm, mainly located in Beilonggang; the maximum horizontal displacement of 1.5 mm is located in the vicinity of Guomeng Town and Longgang Town.
Under the conditions of building loads and current groundwater exploitation, the predicted cumulative maximum ground rebound from 1 January 2021 to 31 December 2030 in Yancheng City is 81.5 mm, and the maximum ground rebound rate is 8.15 mm/a, mainly located in the towns of Dagang and Bicang. The predicted cumulative maximum ground settlement is 3 mm, and the maximum ground settlement rate is 0.3 mm/a, mainly located in the vicinity of Chengbei District and Longgang Town.
Vertically, the submersible aquifer and the group I weak aquifer are the main layers in which subsidence occurs, accounting for 45.5 per cent of the total deformation; the group II and III compressive aquifers are the main layers in which soil rebound is triggered by the current groundwater exploitation, accounting for 47.8 per cent of the total deformation. The recovery of the groundwater level is characterised by rapid growth in the early stage and gradual stabilisation in the later stage.

5. Conclusions

  • In this study, based on the free surface solution method and the incremental solution method for elasto-plastic problems, the Fortran 95 language was used to prepare the three-dimensional coupled numerical analysis programme FGS-3D for groundwater mining and ground settlement. In order for the software to be more operable and to provide a better representation of the calculation results, the Groundwater System Model (GSM) was developed using the Microsoft VisualBasic 6.0 language under the Microsoft Windows 11 operating system environment.
  • On the basis of fully investigating and analysing the hydrogeological and engineering geological conditions of Yancheng City, as well as the soil mechanical properties of the Quaternary strata and the distribution characteristics of the loads of high-rise buildings on the ground, a three-dimensional fully coupled finite-element numerical model of the loads of buildings, groundwater seepage, and deformation of the soil in Yancheng City was established. After identification of the model, the simulation predicted the year-by-year trends of soil compression deformation and horizontal deformation from January 2021 to December 2030 in terms of the superimposed effects of high-rise building loads on the ground and groundwater extraction.
  • Under the action of construction loads and the current groundwater exploitation, the water level of the group III pressurised aquifer in Yancheng City rose by 1.26 m on average in 10 years. In the main groundwater extraction area of the group III pressurised aquifer, two smaller landfall funnels are formed, and the lowest water level of these two landfall funnels is −15 m. The 10-year cumulative maximum ground rebound was 81.5 mm, with a ground rebound rate of 8.15 mm/a. There are three typical ground rebound funnels in Yancheng City, including the large-scale rebound funnel centred on Dagang and Bicang towns, the subscale rebound funnel centred on Louwang and Xuefu towns, and the very small rebound funnel in the city centre.
  • Under the action of construction loads and current groundwater exploitation, the maximum compression of submersible aquifer in Yancheng City from 1 January 2021 to 31 December 2030 is located in the west part of the city; the maximum compression of the weakly permeable layer of the group I clayey soil is also located in the west part of the city; the maximum rebound of the group I compressible aquifer is located in the east part of Dagang Town; the maximum rebound of the weakly permeable layer of the group II clayey soil is located in Beilonggang; the maximum rebound of the group II compressible aquifer is located in Bencang Town and Dagang Town; the maximum rebound of the group III clayey soil is located in Xinxing Town. The maximum resilience of the group II clayey soil weakly permeable layer is located in Beilonggang; the maximum resilience of the group II compressible aquifer is located in Bicang Town and Dagon Town; the maximum resilience of the group III clayey soil weakly permeable layer is located in the north part of the city and Longgang Town; the maximum resilience of the group III compressible aquifer is located in Xinxing Town.

Author Contributions

Methodology, J.D.; software, Z.L.; validation, Y.Z. and C.Z.; formal analysis, J.D.; data curation, J.D.; writing—original draft preparation, J.D.; writing—review and editing, J.D.; project administration, Z.L.; funding acquisition, Y.Z. All authors have read and agreed to the published version of the manuscript.
All authors have read and agreed to the publication of the manuscript.

Funding

This research was funded by Key Laboratory of Earth Fissures Geological Disaster, Ministry of Natural Resources, China (grant number EFGD20210507).

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Gong, S.L. Summary of land subsidence research in Shanghai. Shanghai Geol. 2006, 4, 25–29. [Google Scholar]
  2. Fu, X.G. Cone of Groundwater Depression and Its Development Tendency at Hengshui. Water Resour. Prot. 2001, 1, 24–25. [Google Scholar]
  3. Su, C.; Cui, Y.L.; Shao, J.L. Advances in numerical models of land subsidence based on comprehensive theories. Hydrogeol. Eng. Geol. 2014, 41, 147–152. [Google Scholar]
  4. Yan, T.Z. New mode of forecasting land subsidence regularity. Earth Sci. J. China Univ. Geoscicnces 1989, 14, 181–188. [Google Scholar]
  5. Ma, F.; Yang, F.J. Effect of groundwater exploitation on land subsidence in Tianjin using multiple regression analysis method. Chin. J. Geol. Hazard. Control 2008, 19, 63–66. [Google Scholar]
  6. Zhang, X.L. Analysis and grey prediction of land subsidence variations in Shanghai. Shanghai Geol. 1991, 1, 42–46. [Google Scholar]
  7. Wang, W.; Lu, Y.; Dong, K.G. Grey-Markov Predicting Model to Land Subsidence Evolvement. City Geol. 2008, 3, 41–44. [Google Scholar]
  8. Okumura, T. Analysis of land subsidence in Niigata. In Proceedings of the Land Subsidence Proceeding of the Tokyo Symposium, Tokyo, Japan, 17–22 September 1969; IAHS Publication: Wallingford, UK, 1969; Volume 88, pp. 128–143. [Google Scholar]
  9. R Courant. Variational Method for the Solution of Equilibrium and Vibrations. Bull. Am. Math. Soc. 1943, 49. [Google Scholar]
  10. Ren, R. The preliminary study of relationship between groundwater with drawal and vground subsidence in Cang Zhou city. J. Geol. Hazard. Control 1991, 2, 90–96. [Google Scholar]
  11. Mei, S.B. Discussion on prediction method of surface subsidence in Cang Zhou city. J. Geol. Hazards Environ. Preserv. 1995, 6, 35–43. [Google Scholar]
  12. Lu, L.; Wu, J.C. Bayesian analysis of uncertainties in groundwater numerical simulation. J. Hydraul. Eng. 2010, 41, 264–271. [Google Scholar]
  13. Zhang, X.M. Numerical simulation of gas-waler leakage flow in a two layered coal bed system. J. Hydrodyn. 2009, 21, 692–698. [Google Scholar] [CrossRef]
  14. Zhang, Y.J.; Zhang, Y.J. The numerical analysis of elastic visco-plastic biot’s consolidation to marine soft soil. J. Jilin Univ. (Earth Sci. Ed.) 2003, 33, 71–75. [Google Scholar]
  15. Francesco, R.D. Exact solution of terzaghi’s consolidation equation and extension to two/three-dimensional cases. Appl. Math. 2013, 4, 713–717. [Google Scholar] [CrossRef]
  16. Rogers, J.D.; Chung, J.-W. Applying terzaghi’s method of slope characterization to the recognition of holocene land slippage. Geomorphology 2016, 265, 34–44. [Google Scholar] [CrossRef]
  17. Xu, Y.S.; Shen, S.L.; Cai, Z.Y.; Zhou, G.Y. The state of land subsidence and prediction approaches due to groundwater withdrawal in China. Nat. Hazards 2008, 45, 123–135. [Google Scholar] [CrossRef]
  18. Zhang, Y.; Yu, J.; Gong, X.; Wu, J.; Wang, Z. Pumping-induced stress and strain in aquifer systems in Wuxi, China. Hydrogeol. J. 2018, 26, 771–787. [Google Scholar] [CrossRef]
  19. Wang, Y.; Zhang, F.; Liu, F. Thermo-hydro-mechanical (THM) coupled simulation of the land subsidence due to aquifer thermal energy storage (ATES) system in soft soils. J. Rock Mech. Geotech. Eng. 2024, 16, 1952–1966. [Google Scholar] [CrossRef]
  20. Borja, R.I. Free Boundary Fluid Flow and Seepage Force in Excavation. J. Geo-Tech. Eng. 1992, 1, 125–146. [Google Scholar] [CrossRef]
  21. Zhang, Y.T.; Zhang, W.G. The Boundary Element Method for Seepage Flow in Semi-Infinite Region. J. Hydrodyn. 1981, 4, 8–17. [Google Scholar]
  22. Zienkiewicz, O.C.; Cheung, Y.K. Finite Element in Solution of Field Problems. Engineering 1965, 220, 5710. [Google Scholar]
  23. Ma, Q.S.; Luo, Z.J. Numerical simulation of groundwater exploitation and land subsidence in Cangzhou City. Water Resour. Prot. 2015, 31, 20–26. [Google Scholar]
  24. Luo, Z.J.; Ning, D.; Du, J.J.; Lu, W. Influence of Building Load and Groundwater Exploitation on Land Subsidence in Shengze, Wujiang. J. Jilin Univ. (Earth Sci. Ed.) 2019, 49, 514–525. [Google Scholar]
  25. Luo, Y.; Ye, S.J.; Wu, J.C.; Yang, T.L. Development application of visualization system for numerical simulation of regional land subsidence. J. Huazhong Univ. Sci. Tech. (Nat. Sci. Ed.) 2018, 46, 28–33. [Google Scholar]
  26. Luo, Z.J.; Zeng, F. Finite element numerical simulation of land subsidence and groundwater exploitation based on visco-elastic-plastic Biot’s consolidation theory. J. Hydrodyn. 2011, 23, 615–624. [Google Scholar] [CrossRef]
  27. Vianoa, C.C. The Principle of Rheology in Soil Mechanics. Translated by Du Yupei; Science Press: Beijing, China, 1987. [Google Scholar]
  28. Chen, X.P.; Bai, S.W. The numerical analysis of visco-el astic-plastic Biot ‘s consolidation for soft clay foundation. Chin. J. Geotech. Eng. 2001, 23, 481–484. [Google Scholar]
  29. Luo, Z.J.; Zheng, F.; Tong, X.P. Three-dimensional coupled models of groundwaterseepage and land-subsidence and development of computer program. Resour. Surv. Environ. 2009, 30. [Google Scholar]
  30. Leake, S.A.; Prudic, D.E. Documentation of a Computer Program to Simulate Aquifer-System Compaction Using the Modular Finite-Difference Ground-Water Flow Model; U.S. Geological Survey Techniques of Water-Resources Investigations: Washington, DC, USA, 1991. [Google Scholar]
  31. Fan, X.T.; Ji, G.M. Preconditioned matrix and its stucture technique. J. Cengdu Univ. Technol. (Sci. Technol. Ed.) 2003, 30, 432–435. [Google Scholar]
  32. Luo, Z.J.; Wu, Y.X. Three-dimensional numerical model for groundwater resource evaluation and planin Yancheng. Water Resour. Prot. 2005, 21, 37–41. [Google Scholar]
  33. Jin, W.; Luo, Z.; Wu, X. Sensitivity analysis of related parameters in simulation of land subsidence and ground fissures caused by groundwater exploitation. Bull. Eng. Geol. Environ. 2016, 75, 1143–1156. [Google Scholar] [CrossRef]
Figure 1. Block diagram of the structure of the free surface iteration algorithm.
Figure 1. Block diagram of the structure of the free surface iteration algorithm.
Water 16 02487 g001
Figure 2. Block diagram of the iterative algorithm for viscoplastic-elastic-plastic stress analysis.
Figure 2. Block diagram of the iterative algorithm for viscoplastic-elastic-plastic stress analysis.
Water 16 02487 g002
Figure 3. Block diagram of the fully coupled analysis procedure.
Figure 3. Block diagram of the fully coupled analysis procedure.
Water 16 02487 g003
Figure 4. Functions of software modules.
Figure 4. Functions of software modules.
Water 16 02487 g004
Figure 5. Grid subdivision stereogram of the study area.
Figure 5. Grid subdivision stereogram of the study area.
Water 16 02487 g005
Figure 6. Parameter zoning diagram of the third confined aquifer.
Figure 6. Parameter zoning diagram of the third confined aquifer.
Water 16 02487 g006
Figure 7. Calculated flow field diagram of the third confined aquifer on 31 December 2016.
Figure 7. Calculated flow field diagram of the third confined aquifer on 31 December 2016.
Water 16 02487 g007
Figure 8. Calculated flow field diagram of the third confined aquifer on 31 December 2018.
Figure 8. Calculated flow field diagram of the third confined aquifer on 31 December 2018.
Water 16 02487 g008
Figure 9. Map of cumulative land subsidence from 2016 to 2018.
Figure 9. Map of cumulative land subsidence from 2016 to 2018.
Water 16 02487 g009
Figure 10. Fitted map of ground subsidence at level observation points from 2016 to 2018.
Figure 10. Fitted map of ground subsidence at level observation points from 2016 to 2018.
Water 16 02487 g010
Figure 11. Isocontour map of forecasted water levels of the third confined aquifer on 31 December 2030 (m).
Figure 11. Isocontour map of forecasted water levels of the third confined aquifer on 31 December 2030 (m).
Water 16 02487 g011
Figure 12. Isocontour map of the compressive capacity prediction of the third confined aquifer on 31 December 2030 (mm).
Figure 12. Isocontour map of the compressive capacity prediction of the third confined aquifer on 31 December 2030 (mm).
Water 16 02487 g012
Figure 13. Projected horizontal displacement of the third confined aquifer on 31 December 2030 (m).
Figure 13. Projected horizontal displacement of the third confined aquifer on 31 December 2030 (m).
Water 16 02487 g013
Figure 14. Contour map of the predicted cumulative land subsidence from 2021 to 2030 (mm).
Figure 14. Contour map of the predicted cumulative land subsidence from 2021 to 2030 (mm).
Water 16 02487 g014
Table 1. Comparison table of aquifer group division.
Table 1. Comparison table of aquifer group division.
AquiferAgeRock StratumAquifer
Thicknesses (m)Top Plate Burial Depth (m)Substrate Burial Depth(m)
Submerged and slightly pressurised aquifer groupsHolocene Silt tip group5~40/20~40
Group I pressurised aquiferLate PleistoceneGunnan Group5~2020~4545~70
Group II pressurised aquiferMiddle PleistoceneXiaoxiangzhuang Group20~5050~100130~170
Group III pressurised aquiferEarly PleistoceneFive Teams Township Group10~60150~180190~270
Table 2. Summary of parameter values for each parameter partition in the group III pressurised aquifer.
Table 2. Summary of parameter values for each parameter partition in the group III pressurised aquifer.
PartitionsCoefficient of Permeability in the Direction of the Main Axis (m·d−1)Deformation Modulus (MPa)Poisson’s RatioCohesion (KPa)Friction Angle (°)Expansion Angle (°)Severity (KN·m−3)Effective Porosity
KxxKyyKzzEνcφψΓn
4413131.3450.395.834.8020.180.435
4512121.2460.416.234.2020.540.432
4618181.8500.4633.7021.620.426
4720202480.385.933.4020.350.42
4816161.6460.366.532.8021.240.433
4912121.2530.396.432.5021.330.428
5014141.4480.376.133.4021.170.419
51990.9490.415.835.5020.930.413
528.58.50.85510.435.736.1020.540.417
53330.3470.455.337020.360.437
544.54.50.45450.46537.5020.220.422
Table 3. Table of inflow and outflow from the group III confined aquifer.
Table 3. Table of inflow and outflow from the group III confined aquifer.
VintagesInflow (m3·d−1)Outflow (m3·d−1)
2016112,971.3107,914.38
2017123,832.22117,363.11
2018114,567.3590,794.41
2019109,878.8188,460.83
202091,473.8870,632.2
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Du, J.; Zhang, Y.; Luo, Z.; Zhang, C. Prediction of Ground Subsidence Induced by Groundwater Mining Using Three-Dimensional Variable-Parameter Fully Coupled Simulation. Water 2024, 16, 2487. https://doi.org/10.3390/w16172487

AMA Style

Du J, Zhang Y, Luo Z, Zhang C. Prediction of Ground Subsidence Induced by Groundwater Mining Using Three-Dimensional Variable-Parameter Fully Coupled Simulation. Water. 2024; 16(17):2487. https://doi.org/10.3390/w16172487

Chicago/Turabian Style

Du, Jingjing, Yan Zhang, Zujiang Luo, and Chenghang Zhang. 2024. "Prediction of Ground Subsidence Induced by Groundwater Mining Using Three-Dimensional Variable-Parameter Fully Coupled Simulation" Water 16, no. 17: 2487. https://doi.org/10.3390/w16172487

APA Style

Du, J., Zhang, Y., Luo, Z., & Zhang, C. (2024). Prediction of Ground Subsidence Induced by Groundwater Mining Using Three-Dimensional Variable-Parameter Fully Coupled Simulation. Water, 16(17), 2487. https://doi.org/10.3390/w16172487

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop