Land subsidence characteristics and disaster prevention in the Tongzhou area, Beijing

Accurately assessing the impact of groundwater exploitation on land subsidence can provide scientific support for decision-makers. This article establishes a groundwater–subsidence model in a typical area using the land subsidence characteristics of Tongzhou, China, and subsequently classifies a land subsidence early warning zone based on the land subsidence rate. Results show that land subsidence occurred throughout the whole Tongzhou District. The land subsidence that developed to the west of Liuzhuang, Tongzhou city and Taihu was the most serious The maximum annual subsidence rate reached 120 mm a⁻¹ under the current groundwater extraction conditions, and the early warning level for land subsidence reached its highest level. If the pumping of groundwater was reduced by 50 % in the second and the fourth aquifers and by 60 % in the third aquifer, the land subsidence early warning level would be largely reduced and would meet the requirements for land subsidence control.

1 Introduction

Land subsidence, which is a slowly developing geological disaster, has become a global geological problem characterized by a long duration, slow generation, a wide influence range, a complex mechanism, and great difficulty with respect to prevention and control (Martinez-Noguez and Hinkelmann, 2015) There are currently more than 60 countries and regions in the world suffering from land subsidence (Yuan et al., 2012). UNESCO, IAHS and other groups have carried out a series of international research projects on land subsidence with research topics including monitoring land subsidence, the mechanism behind land subsidence and its prevention. In Beijing, subsidence mainly develops in Chaoyang, Haidian, Shunyi, Tongzhou, Changping and Daxing. The historical maximum land subsidence rate has reached 159.6 mm a⁻¹, and the land subsidence rate is currently still rapidly developing (Zhao et al., 2018a).

The effect of land subsidence on urban planning, development and economic construction has become more prominent. Land subsidence has induced ground fissures in complex geological environments, including Shunyi and Changping, which have subsequently caused damage to infrastructure and reduced the value of land resources (Zhao et al., 2018b). Land subsidence reduces the ground elevation, causing the lifting, tilting and cracking of deep-water well pipes. It has also threatened the operation of rail transit (Liu et al., 2016).

This paper aims to establish an early warning mechanism with respect to land subsidence for disaster prediction and to reduce and control the multiple disasters induced by this phenomenon. This work utilized Tongzhou as the study area. A regional soil–water coupling model based on a three-dimensional groundwater model and a one-dimensional vertical soil deformation model was established. Based on the simulation outcomes from this coupled model, the warning level for land subsidence and possible prevention and control measures were then provided. This study will supply scientific support to ensure the safety of major infrastructure projects, buildings and urban operations.

2 Characteristics of aquifer structure and land subsidence

The study area is located in the Tongzhou and Chaoyang districts, which encompass a total area of about 300 km². Land subsidence in this area is serious.

2.1 Aquifer structure

The sediments of the Quaternary strata in Tongzhou are formed from the Yongding River alluvial fan and the Chaobai River alluvial fan. The Quaternary strata is less than 300 m thick in the Maji Bridge, Taihu and Yongshun regions, and the thickness increases gradually to the east and southeast, reaching more than 500 m. The Quaternary strain shows gray, gray black and yellow colors. The lithology is clay, clay silt, sandy silt, fine sand, medium coarse sand or gravel. The alternating distribution of sand and clay layers in the strata form the multi-structure.

The study area encompasses the Yongding River, Chaobai River and Wenyu River alluvial plain. Within a depth of 30 m below the ground surface, the soil type sequence comprises an artificial soil layer, a Quaternary new sedimentary layer and a general Quaternary alluvial deposit. There are also two thick sandy soil layers; the sandy soil content increases from west to east, and the depth of the sandy soil layers became shallower.

There is an aquitard at a depth of about 100 m based on an analysis of borehole data. The Quaternary stratum above 100 m is defined as shallow aquifer, whereas that below 100 m is defined as confined aquifer. The lithology is composed of alternative sand and clay, and hydraulic connections exist between the two abovementioned aquifers.

The shallow aquifer can be divided into a phreatic aquifer and a shallow confined aquifer; the first aquifer system has a bottom depth of about 40–50 m, and the second aquifer system has a bottom depth of approximately 80–100 m. The aquifer is mainly composed of fine sand, followed by medium sand and coarse sand, and gravel in local areas. The aquifer generally has five to eight layers with a cumulative thickness ranging from 30 to 50 m. The aquifer recharge conditions are good in the central and northern parts of Tongzhou.

2.2 Subsidence characteristics

2.2.1 Characteristics of horizontal subsidence

The subsidence boundary for Tongzhou has gradually merged with the boundary of subsidence in the Chaoyang District, including the Jinzhan and Heizhuanhu regions (Luo, 2019). Currently, the land subsidence in Tongzhou is still rapidly developing, and the subsidence center is bordered by the downtown area of Tongzhou, Zhangjiawan, Taihu and Ciqu. The maximum cumulative subsidence and subsidence rate in the subsidence center until 2016 was about 1700 mm and 110 mm a⁻¹, respectively. Cumulative subsidence of approximately 500 mm has been observed in the areas west of Songzhuang in Tongzhou and west of Taihu, whereas a higher cumulative subsidence of 1000 mm has been found to the west of Liuzhuang, Tongzhou and Taihu. In Tongzhou and at the junction of Taihu and the Chaoyang, a cumulative subsidence as high as 1500 mm has been recorded.

2.2.2 Characteristics of vertical subsidence

The Zhangjiawan extensometer station began operations on 26 June 2008, and it recorded a cumulative subsidence of 139.128 mm from June 2008 to 26 December 2015. The time series subsidence curve at each monitoring layer showed an increasing trend (Luo, 2019). Land subsidence in the deep stratum has been higher than that in the shallow stratum since 2008.

The contribution of different layers to the cumulative land subsidence showed that the stratum at a depth below 126 m accounted for about 93.4 % of the land subsidence in the Zhangjiawang region. The deformation of the layers at depths greater than 257 m accounted for 54.01 % of the total compression value. The deformation of the strata at a depth between 126 and 193 m and between 193 and 257 m accounted for 27.70 % and 11.69 % of the total compression value, respectively, and the strata at depths of 49–90 m, 30–49 m and 2–30 m accounted for 1.99 %, 6.30 % and 3.18 %, respectively. Notably, the stratum at a depth of 90–126 m was uplifted (Luo, 2019).

3 Groundwater–subsidence numerical model

3.1 Conceptual model

The Quaternary strata in the study area can be divided into the Holocene first aquifer, the Upper Pleistocene second aquifer, the Middle Pleistocene third aquifer and the Lower Pleistocene fourth aquifer. This study combined the aquifer and aquitard of Quaternary strata and regarded the clay layer lens and aquifer as a whole groundwater system. The direction of groundwater flow is mainly horizontal. The groundwater flow pattern is set as a three-dimensional unstable flow.

The model boundary was set as the flow rate condition on the base of the historical groundwater flow field. The source/sink terms included river channel infiltration, field irrigation infiltration, atmospheric precipitation, evaporation and pumping. The depth of the model bottom was about 350 m, which contained all of the groundwater exploitation layers.

Table 1 Classification of land subsidence warning levels.

Table 2 The different reduction scenarios of groundwater exploitation.

3.2 Groundwater flow numerical model

The classic three-dimensional saturated flow equation is given in Eq. (1).

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{ll}S\frac{\partial h}{\partial x}=\frac{\partial \left(K_{\mathrm{cp}}\frac{\partial h}{\partial x}\right)}{\partial x}+\frac{\partial \left(K_{\mathrm{cp}}\frac{\partial h}{\partial x}\right)}{\partial y}+\frac{\partial \left(K_{\mathit{\nu }}\frac{\partial h}{\partial x}\right)}{\partial z}+\mathit{\epsilon }& x,y,z\mathit{ϵ}\mathrm{\Omega },t\ge \mathrm{0}\\ \mathit{\mu }\frac{\partial h}{\partial t}=K_{\mathrm{cp}}{\left(\frac{\partial h}{\partial x}\right)}^{\mathrm{2}}+K_{\mathrm{cp}}{\left(\frac{\partial h}{\partial y}\right)}^{\mathrm{2}}& \\ +K_{\mathit{\nu }}{\left(\frac{\partial h}{\partial x}\right)}^{\mathrm{2}}+\frac{\partial h}{\partial x}\left(K_{\mathit{\nu }}+p\right)+p& x,y,z\mathit{ϵ}{\mathrm{\Gamma }}_{\mathrm{0}},t\ge \mathrm{0}\\ h\left(x,y,z,t\right)\left|{}_{t=\mathrm{0}}\right.=h_{\mathrm{0}}& x,y,z\mathit{ϵ}\mathrm{\Omega }\\ \left.h\right|{\mathrm{\Gamma }}_{\mathrm{1}}=h_{\mathrm{1}}& x,y,z\mathit{ϵ}{\mathrm{\Gamma }}_{\mathrm{1}},t\ge \mathrm{0}\\ K_n\frac{\partial h}{\partial n}\left|{\mathrm{\Gamma }}_{\mathrm{2}}\right.=q\left(x,y,z,t\right)& x,y,z\mathit{ϵ}{\mathrm{\Gamma }}_{\mathrm{2}},t\ge \mathrm{0}\\ \frac{h_r-h}{\mathit{\sigma }}-K_n\frac{\partial h}{\partial n}\left|{\mathrm{\Gamma }}_{\mathrm{3}}\right.=\mathrm{0}& x,y,z\mathit{ϵ}{\mathrm{\Gamma }}_{\mathrm{3}},t\ge \mathrm{0}\\ \frac{\partial h}{\partial n}\left|{\mathrm{\Gamma }}_{\mathrm{4}}\right.=\mathrm{0}& x,y,z\mathit{ϵ}{\mathrm{\Gamma }}_{\mathrm{4}},t\ge \mathrm{0}\end{array}\right\},\end{array}$$

where Ω is the simulation domain; h is the hydraulic head (m); K_cp, K_v and K_n represent the hydraulic conductivity tensor with the principal components (m d⁻¹); S is the specific storage; μ is the specific yield; σ is the resistance coefficient of the river bed, $\mathit{\sigma }=L/K_s$, where L and K_s are the thickness and vertical permeability coefficient (m d⁻¹) of the river bed, respectively; ε is the source/sink term (d⁻¹); p is evaporation and precipitation (d⁻¹); h0 is the initial water level (m); Γ0 is the upper boundary of the domain; Γ1 is the Dirichlet boundary; Γ2 is the flow rate boundary; Γ3 is the mixed boundary condition; Γ4 is the bottom boundary of the domain; $\stackrel{\mathrm{̃}}{n}$ is the normal direction of the boundary; and $q(x,y,z,t)$ is the per unit area of the flow rate (m d⁻¹), where inflow is positive and outflow is negative.

3.3 Soil displacement model

This paper assumes that the displacement occurred in the vertical direction, considering the relative small deformation in the horizontal direction compared with that in the vertical direction. The change in the groundwater level enabled the calculation of the amount of the soil deformation as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& -\mathrm{\Delta }h{S}_{\mathrm{sk}}=\mathrm{\Delta }b\\ & S_{\mathrm{skv}}={\displaystyle \frac{\mathrm{0.434}{C}_{\mathrm{c}}{\mathit{\gamma }}_{\mathrm{w}}}{{\mathit{\sigma }}^{\prime }(\mathrm{1}+e_{\mathrm{0}})}}({\mathit{\sigma }}^{\prime }>Pc)\\ & S_{\mathrm{ske}}={\displaystyle \frac{\mathrm{0.434}{C}_{\mathrm{r}}{\mathit{\gamma }}_{\mathrm{w}}}{{\mathit{\sigma }}^{\prime }(+e_{\mathrm{0}})}}({\mathit{\sigma }}^{\prime }<Pc),\end{array}\end{array}$$

where Δh is the change in the groundwater level; S_sk is the skeletal specific storage of the aquitard (S_ske is the skeletal elastic storage, and S_skv is the skeletal inelastic storage); Δb is the change in the thickness of the aquitard; C_c is the compression index; C_r is the expansion index; γ_w is the bulk density of water; σ^′ is the effective stress; and e₀ is the initial porosity ratio.

Once the deformation of each aquifer was calculated, the cumulative value of displacement could be summarized. The flow model and the soil deformation model can be coupled via the flow term:

$$\begin{array}{}\text{(3)}& q_i^m={\displaystyle \frac{S_k^m}{\mathrm{\Delta }{t}^m}}(h^m-H^{m-\mathrm{1}})+{\displaystyle \frac{S_{\mathrm{ke}}}{\mathrm{\Delta }{t}^m}}(H^{m-\mathrm{1}}-h^{m-\mathrm{1}}),\end{array}$$

where m denotes the mth time step, i denotes the ith grid, q is the flow of the unit per unit time, Δt^m is the length of the time step m, H^m−1 is the pre-consolidation water level at the end of the time step m−1, h is the calculated water level and $S_k^m$ is determined by

$$\begin{array}{}\text{(4)}& {\displaystyle }{S}_k^m=\left\{\begin{array}{c}{S}_{\mathrm{ke}},h^m>H^{m-\mathrm{1}}\\ S_{\mathrm{kv}},h^m\le H^{m-\mathrm{1}}\end{array}\right..\end{array}$$

3.4 Model calibration

The model was vertically divided into four layers and horizontally divided into 2500 grids (50 in the x direction and 50 in the y direction) with a grid size of 400 m×400 m. The simulation time was from January 2005 to December 2015. There were a total of 131 stress periods with an interval of 1 month. Each period included several time steps, and the error of each iteration was strictly controlled.

The trial–correction method was adopted to identify the hydrogeological parameters, boundary values and other equilibrium terms by fitting the historical groundwater level with the long-term monitored data from the same period. During the process, the simulated volumes of recharge and discharge were also compared with the monitored data. The monitored flow field and simulated flow field in December 2015 were compared for validation. The simulated flow field generally reflects the groundwater level distribution. The groundwater level of the long-term monitored data was fitted with the simulated outcomes, which reflected that the two curves were approximately consistent.

4 Early warning and prevention of land subsidence

The early warning areas and warning levels with respect to land subsidence were differentiated (the latter into three levels) by referring to the “Beijing Land Subsidence Prevention and Control Plan (2013–2020)”, the “Beijing–Tianjin–Hebei Land Subsidence Integrated Prevention and Control Master Plan” and the “Surface Subsidence Investigation and Monitoring Code”. Land subsidence was effectively monitored by InSAR, leveling and GPS. The relations between the subsidence rate and the early warning threshold were determined based on these real-time monitoring data. If the subsidence rate was larger than the threshold, the warning level was given.

In this study, we modeled the subsidence rate in 2020 using the unchanged groundwater exploitation conditions from 2015. The groundwater–subsidence coupling model was utilized with identified and corrected parameters. The results show that the original land subsidence center rapidly developed, and the maximum annual subsidence rate reached 120 mm a⁻¹.

Three kinds of scenarios regarding reductions in the groundwater exploitation were simulated in 2016 to access which scenarios would be effective to reduce the land subsidence from 2016 to 2020 (Table 2). Based on the model predictions, the subsidence rate in scenario 1 is greater than 15 mm a⁻¹, and that in scenario 2 is generally less than 15 mm a⁻¹. However, the subsidence rate for Zhangjiawan, Taihu and Heizhuanghu (the subsidence center) is still greater than 30 mm a⁻¹. Based on the results of scenario 3, the second and fourth aquifer exploitation was set to be reduced by 50 %, and the third aquifer was reduced by 70 % in Zhangjiawan, Taihu, and Heizhuanghu. The modeling results show that the maximum annual subsidence rate is about 30 mm a⁻¹, and the regional subsidence rate is less than 15 mm a⁻¹, which is lower than the warning threshold in 2020.

5 Results

The three main results of this study are as follows:

1. The Tongzhou district has developed serious land subsidence. The subsidence of the strata at depths below 126 m has dominantly contributed to the total subsidence. Land subsidence in the area comprising Liuzhuang, Tongzhou and downtown Taihu bordering the Chaoyang District is the most serious. The shallow strata show little subsidence, and the stratum from a depth of 90 to 126 m was even uplifted.

2. We established the groundwater flow coupled to a soil deformation model, which can be used to simulate the dynamic features of the groundwater level and land subsidence. Based on the model, land subsidence early warning zones and level were given. The model predicted that the regional and central subsidence rates would reach the highest warning level by 2020 under the current groundwater exploitation conditions.

3. Using alternative projects, including rainwater, reclaimed water, south–north water transfer and homemade well replacement, the decrease in the groundwater level could be restricted and land subsidence would be relieved. The land subsidence warning level would be below the minimum threshold warning level if the pumping volume from the second and fourth aquifers was reduced by 60 % and the pumping volume from the third aquifer was reduced by 70 % the Zhangjiawan, Taihu and Heizhuanghu region, and the pumping volume from the second and fourth aquifers was reduced by 50 % and pumping from the the third aquifer was reduced by 60 % in other regions with respect to the current groundwater exploitation conditions.