Integrating SBAS-InSAR and Random Forest for Identifying and Controlling Land Subsidence and Uplift in a Multi-Layered Porous System of North China Plain

Abstract

Controlling groundwater table decline could mitigate land subsidence and induced environmental hazards in over-explored areas. Nevertheless, this becomes a challenge in the multi-layered porous system as (in)elastic deformation simultaneously occurs due to vast spatiotemporal variability in the groundwater table. In this study, SBAS-InSAR was used to estimate annual land deformation during 2017–2022 in a specific region of North China Plain, in which aquifers are composed of many layers of fine-grained compressible sediments and the groundwater table has experienced a prolonged decline. The random forest (RF) was applied to establish the nonlinear relationship between accumulated deformation and its potential driving factors, including the depth to the groundwater table (GWD) and its change rate, and the compressible sediment thickness. Results show that the marked subsidence and uplift co-exist in the region even though the groundwater table has risen widely since the South–North Water Diversion Project. The land subsidence is attributed to inelastic compaction of the thick compressible deposits in depression cone centers, where the GWD is over 40 m and 90 m in the shallow and deep aquifers, respectively. In contrast, the marked uplift is primarily attributed to fast rising of the groundwater table (e.g., −2.44 m/a). The RF predictions suggest that, to control the subsidence, the GWD should be less than 20 and 70 m in the shallow and deep aquifers, respectively, and the rising rate of the GWD should increase to 2–5 times of current rates in the depression cones. To mitigate the marked uplift, the rising rate of the GWD should reduce to 1/2–1/5 of the current rates in the shallow aquifers. The uneven deformations of sediments in the depression cone centers and uplift in their boundaries may exacerbate geohazards. Therefore, it is vital to implement appropriate governance of groundwater recovery in the multi-layered porous system.

1. Introduction

Land deformation, including subsidence and uplift, could result in serious problems, such as increasing the risk of flooding in coastal areas, damaging buildings and infrastructures, destructing local groundwater systems, generating tension cracks on land, and reactivating faults [1,2]. Land subsidence has typically resulted from over-exploitation of groundwater and underground space development [3,4]. Roughly 73% of the mapped subsidence occurs over cropland and urban areas due to subsurface water overdraft [5]. Reducing and ceasing groundwater exploration or artificial recharge could mitigate land subsidence and even lead to ground rebound or uplift.

Even though the mechanism by which rock and sediments deform and compact under the influence of a change in pore water pressure is well understood, land deformation in a multi-layered porous aquifer system is complicated by the combined influences of the effective vertical and horizontal stresses over several layered deposits. In particular, a long-term decline in the groundwater table produces a disturbance that propagates its effect in space and time through the geological medium [6]. Moreover, compared to land subsidence, land uplift is a much less observed and recognized event [7]. The fundamental understanding of land uplift at present is poor. Therefore, accurately predicting the subsidence and uplift rates and identifying the dominant influencing factors in the multi-layered porous system are vital for controlling groundwater explorations and preventing induced geohazards.

Effective monitoring is an initial step to estimate the extent of land deformation. Land deformation is conventionally monitored by geodetic techniques, including the Global Navigation Satellite System (GNSS), bedrock and layering mark surveying, and leveling measurements, although these techniques are limited by precise measurements and a small monitoring range. Satellite technologies, such as Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR), allow regional-scale mapping and monitoring of land deformation with lower costs [8]. In particular, the monitoring of large spatial extent and temporal development of land deformation is possible with the application of time series radar interferometry by analyzing a series of stacked radar interferometric data with continuous temporal baselines [9,10,11,12]. For instance, Galve et al. [13] performed SBAS-INSAR calculations based on the European Space Agency’s Geohazard Exploitation Platform (GEP) on 25 ENVISAT satellite images of southern Spain from 2003 to 2008, and the ground instabilities were detected. Du et al. [14] effectively identified the surface subsidence and uplift in Eastern Tianshan Mountain utilizing SBAS-InSAR from 25 December 2017 to 2 January 2021. Furthermore, Nayak et al. [15] extracted ground displacement in Morocco for the year 2023 based on SBAS to elucidate the impact of earthquakes on the ground.

The ground subsidence magnitude, time of occurrence, and extent of the area depend on a large number of factors including the amount of groundwater withdrawn, the pore pressure decline, the depth and permeability of the pumped formation, the geomechanical properties of the aquifer, and the information of overburden [7]. Over-exploitation of groundwater causes subsidence because it increases effective stress in the aquifer system and results in vertical compaction [16]. If the groundwater head declines below the pre-consolidation head, producing effective stresses greater than the maximum historical stress, inelastic deformation occurs and the sediments will undergo irreversible compaction. The inelastic (non-recoverable) compaction is largely concentrated in the aquitards (compressible silt and clay deposits) [17,18]. The uplift can only arise as a result of water pressure rebound in elastically compaction areas [19]. Therefore, whether the cessation of groundwater over-exploitation could lead to a ground rebound depends on the extent of the historical water level and the composition of sediments in the aquifer system.

The theoretical process-based models can be used to predict subsidence and uplift rates. However, the models require extensive geomechanical and hydrogeologic datasets, and knowledge of the temporal evolution of head changes driving nonlinear subsidence processes, which hamper applications in a large-scale multi-layered aquifer system. In particular, as the groundwater table changes from a decline to a rise in the multi-layered aquifer system, the shift from draining to re-wetting of pore space results in a change in the geostatic stress on the underlying sediments. Because of the high-dimensional and nonlinear problem of variables, identifying the effects of multiple variables on land deformation is also difficult to achieve by traditional statistics methods [4].

In recent years, machine learning (ML), with its advantage of dealing with multidimensional data, has been frequently used for modeling systems with nonlinear correlation structures between the dependent and independent variables [20]. In particular, the random forest (RF) algorithm, as an ensemble machine learning method based on a popular tree-based ensemble learning algorithm, obtains excellent accuracy when handling multivariate complex nonlinear problems [5]. RF has been widely used to model complex nonlinear relationships, including the problems with aquifers. For instance, Chen et al. [21] utilized RF to explore the relationship between land subsidence and groundwater in different aquifers in Beijing Plain during 2011–2018, which indicated groundwater in the second confined aquifer had the biggest impact on the subsidence of all aquifers. Fu et al. [22] applied RF to investigate the spatial distribution of temporal correlation, and it could explain 72% and 60% of the spatial variance between the ground deformation and critical head/groundwater level, respectively.

The North China Plain (NCP) is one of the regions typified by the long-time evolutionary process of land subsidence [1]. The NCP is a large-scale down-faulted rift basin formed in the late Paleogene and Neogene. It was then modified by alluvial and river deposits, resulting in multi-layered aquifer systems composed of numerous lenses or layers of fine-grained compressible sediments. The NCP has suffered from prolonged groundwater storage depletion and land subsidence since the 1970s due to the rapidly growing demand for agricultural, industrial, and domestic water uses [23]. The maximum decline of the groundwater table even exceeded 100 m in the center of the groundwater depression zone in the NCP [24]. In order to reduce over-drafting of the aquifer system and mitigate the negative effects of land subsidence, the South-to-North Water Diversion Project began operation to reduce groundwater utilization in 2014. Since then, the groundwater table has risen, and ground subsidence has been reduced and even rebounded in some areas. Ground subsidence has been widely investigated, mostly focusing on estimates of the subsidence rates and groundwater storage in relation to the prolonged groundwater decline [25,26,27]. The ground uplifts in some areas have been reported based on limited observation of borehole extensometer stations [28] and the monitoring by remote sensing of SAR datasets during the early period of the South-to-North Water Diversion Project and the groundwater exploitation restriction policy [23,29]. Because of the time-lag between land deformation and groundwater table variations, these investigations cannot sufficiently explore where the ground subsidence could be reduced or shifted to rebound with the consecutive recovery of the regional groundwater level.

The purpose of this paper is to identify subsidence and uplift areas, and to explore their relations to the driving factors of changes in groundwater table and compressible sediments in a multi-layered porous aquifer system in the Cangzhou–Hengshui region of the NCP. The land deformation after the regional rising of the groundwater table from 2017 to 2022 due to a massive water division was extracted using the SBAS-InSAR technique. The remarkable subsidence and uplift areas were identified based on the deformation rate being larger or smaller than a threshold (e.g., $\pm$10 mm/a). The RF method was then applied to identify the main driving factors of subsidence and uplift rates by training decision trees to establish multi-dimensional relationships between land deformation and its influencing factors. Finally, the trained RF was used to generate possible deformation rates resulting from changes in the driving factors, and to select appropriate combinations of driving factors that could alleviate ground deformation rates. Using the prediction method and results, one can understand the subsidence evolution under human efforts to recover the groundwater level and to make regional water management policies considering the multi-layered porous aquifer system.

Study Area

The North China Plain (NCP) has a total area of about 140,000 km², and comprises the megalopolis of Beijing and Tianjin as well as the Hebei Province of China (Figure 1). The Cangzhou–Hengshui region is located in the east of Hebei Province (37°03′N~38°57′N and 115°10′E~117°52′E), with an area of about 2.3 × 10⁴ km² (Figure 1). This region is one of the most serious subsidence sites in the NCP, where the accumulative subsidence was up to −3 m from 1971 to 2003 [25]. The climate is dominated by the continental monsoon, with annual precipitation ranging from 500 to 600 mm, and 50 to 80% of the total occurring in the summer monsoon months from July to September.

The main regional strata are unconsolidated Quaternary sediments with thickness ranging from 350 m to over 550 m (Figure 2). The aquifer system is composed of sand and gravel interbedded with fine sand and silt layers, which can be divided into four major aquifer groups (I–IV in Figure 2) [30]. In the Cangzhou–Hengshui region, the thickness of the aquifers (I, II, III, and IV) is 20–63 m, 108–194 m, 105–232 m, and 71–279 m, respectively. The lenses or layers of fine-grained clay and silt occupy a high percentage of the entire deposit thickness, e.g., over 50% for the confined aquifer III. The top aquifer is unconfined, while the other three are (semi-)confined. The two upper aquifers (I and II) are hydraulically connected in most areas due to extensive well drillings and the discontinuity of the aquitards. Thus, the two upper aquifers (I and II) are usually merged into an unconfined or a semi-confined aquifer.

The NCP is considered one of China’s productive agricultural regions. It has experienced a massive expansion of irrigated agriculture since the 1960s [31,32]. Because of three marine transgressions in the coastal area of Bohai Bay since the Late Pleistocene, groundwater in the upper aquifer system is saline or brackish in the eastern regions, including the Cangzhou–Hengshui region. Because the bottom depth of the saline/brackish water increases from 40 m in the west to more than 300 m in the east [23], the freshwater utilizations mainly come from the upper aquifers (I and II) in the northwest of the study area and primarily pumped water from the confined aquifer systems (III and IV).

2. Data and Methods

2.1. Synthetic Aperture Radar Data

The Synthetic Aperture Radar (SAR) images were obtained from the Sentinel-1 constellation, which was developed for the Copernicus joint initiative of the European Commission and the European Space Agency. The Sentinel-1 satellite has a revisit period of 12 days. The single-look complex SAR datasets were acquired in interferometric wide swath (IW) mode, with Terrain Observation with Progressive Scan (TOPS) imaging geometry. These datasets are imaged in C-band at an operating frequency of 5.405 GHz in the right-looking geometry. In this study, 83 ascending Sentinel-1 images (Figure 1) were collected from January 2017 to December 2022, with an average time interval of 24 days.

2.2. SBAS-InSAR

SBAS-InSAR based on the Differential Interferometric Synthetic Aperture Radar (DInSAR) technique was originally proposed by Berardno [33]. It was improved by combining the differential interferograms generated from pairs characterized by small orbital baselines to limit spatial correlation. The SAR images in Track-121 and Track-116 (shown in Figure 1) at the same monitoring time were spliced to obtain an image that covers the entire study area. Interferogram generation was completed by sampling and matching the SAR images according to the baseline threshold. Then, the topographic phase was removed from the interferometric phase with the 30 m resolution SRTM Digital Elevation Model (DEM). High-pass spatial filtering and low-pass temporal filtering steps were applied to estimate and remove atmospheric effects. To solve the series of interferograms characterized by small spatiotemporal baselines, the least squares method (LSM) was applied to retrieve deformation information on each small spatiotemporal baseline, and the singular value decomposition (SVD) algorithm was utilized to solve multiple baseline subsets to obtain the line of sight (LOS) displacement in the entire time series. Finally, the results of SBAS-InSAR in each point of the frames were transformed from LOS to the vertical direction by:

$$d_{LOS}=\mathit{cos}{\theta }^i·d_U$$

(1)

where $d_U$ represents surface deformation of the vertical, $d_{LOS}$ represents surface deformation in the LOS direction, and $\mathit{cos}{\theta }^i$ represents the radar incidence angles at the i-th point of the frames.

In this study, the groundwater-induced land deformation is usually dominated by vertical land deformation, and the horizontal deformation could be negligible. The vertical deformation was calculated according to the LOS direction of InSAR-retrieved deformation. The negative values of pixel deformation represent ground subsidence, while the positive values represent an uplift.

In order to provide an appropriate sample dataset capacity for RF modeling, the fine resolution values of land deformation at the pixels (10 m $\times$ 10 m) were averaged in a coarse resolution of grids (500 m $\times$ 500 m). The gridded deformation at intervals of 24 days was accumulated into monthly and annual values from 2017 to 2022.

2.3. The Predictor Variables

The predictor variables can be combined into two groups: (1) the depth to the groundwater table (GWD) and its change rate, which reflects the effects of the pressure head on effective stress; and (2) thickness of the compressible sediments (silt and clay), which quantifies the inelastic (non-recoverable) compaction of the sediments. For the three aquifer groups (I & II, III, and IV), there are nine predictors in total.

In the study region, the groundwater level or depth was observed at 839 monitoring wells, among which 370, 212, and 257 wells were set in the aquifers of I & II, III, and IV, respectively. The recording interval was mostly 10 days (the 1st, 11th, and 21st) in each month from 2017 to 2022. Additionally, there are core-analysis data at 143 boreholes (Figure 1) at depths ranging from 25 to 580 m in the study region. According to the data, the composition of Quaternary sediments can be classified as sand, silt or silty clay, and clay. Since clay and silty clay have relatively high plasticity, they are identified as highly compressible sediments [31].

The spatial distributions of the GWD and the thickness of sand, silt, and clay layers were interpolated using the Kriging algorithm by leveraging the observed water table at the monitoring wells and the measured thickness at the boreholes. The interpolated spatial distributions were then used to extract the gridded values with a spatial resolution of 500 m $\times$ 500 m, and to obtain the thickness and proportion of compressible sediments in each aquifer.

The change rate of a variable (e.g., land deformation or GWD) was obtained from the slope of the fitted linear function for the annual series.

2.4. Random Forest

2.4.1. Construction of Random Forest

The random forest (RF) model was used to analyze the nonlinear relationship between land deformation and its associated factors, i.e.,

$$Y=F\left(x_1, x_2\dots \dots x_n\right)$$

(2)

where $Y$ represents the predicted values of ground deformation (the accumulative deformation in this study), $x_i$ represents the ith variable (i.e., mean depth and annual change rate of the GWD, and compressible deposit thickness in each aquifer), n is the number of varibles, and $F\left(\right)$ represents the high-dimensional and nonlinear structure within the random forest model.

As an ensemble model, RF combines both regression trees and the boosting technique, which allows the production of a large number of simple tree models and then combines them into a final optimized model [34]. The RF is useful for handling different types of predictor variables without prior data transformation or outlier elimination [35]. Random forests for regression are formed by growing a certain number of decision trees based on the CART algorithm [36]. Each regression tree is trained based on the generated dataset of randomly selected variables. The final prediction of an RF predictor is calculated by taking the average over the results of all the decision trees.

2.4.2. Validation of RF Model

We used a ten-fold cross-validation method to validate the RF model. That is, in the process of model implementation, 70% of the samples (8649 cells of the subsidence area and 4396 cells of the uplift area) were randomly selected for RF regression training and the remaining 30% (3708 cells of the subsidence area and 1319 cells of the uplift area) were used for accuracy verification. The ratio of 70% of the dataset for training and 30% for testing is considered the best ratio for training and validating [37,38]. Since the training dataset is randomly sampled each time, we repeated the process ten times and took the best value of accuracy to ensure the reliability of the validation results.

To assess the accuracy of the fitted models, the coefficient of determination (R²) was used:

$$R^2=1-\frac{u}{v}$$

(3)

where $u$ represents the residual (between the predicted value and actual value) sum of squares, and $v$ represents the total sum of squares. $R^2$ ranges from zero to one, where $R^2=1$ means the predicted value perfectly captures the observations.

In addition, mean absolute error (MAE), and root mean square error (RMSE) were also calculated to validate the RF models as follows:

$$\mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{N}\sum _{i=1}^N\left|{Yo}_i-{Yc}_i\right|$$

(4)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{N}\sum _{i=1}^N{({Yo}_i-{Yc}_i)}^2}$$

(5)

where N is the number of samples, and ${Yo}_i$ and ${Yc}_i$ are observed and predicted values at the ith grid, respectively.

2.4.3. Extraction of Importance Scores

The RF algorithm allows assessment of the relative importance of the input variables, or the variable importance. The variable importance can be estimated by the permutation importance as follows [34]:

$$I_x={\sum }_{k=1}^K\left[\frac{1}{K}({MSE}_k^{x_{perm}}-{MSE}_k)\right]$$

(6)

where $I_x$ is the importance of variable $x$, K is the number of trees in the forest,${MSE}_k^{x_{perm}}$ is the estimation error with predictor $x$ being eliminated for the kth decision tree, and ${MSE}_k$ is the forecasting error with all predictors included in the kth decision tree. A higher score of $I_x$ means the variable $x$ is more relevant to the regression. Usually, the importance scores of all the input variables are summed to 100.

2.4.4. Threshold of Influencing Factors for Deformation Control

Based on the prior distribution of influencing factors (as delineated in

Table 1. Characteristics of land deformation and its influencing factors, including the groundwater table depth (GWD) and its change rate, and compressible deposit thickness in the classified areas.

Classified Area		Area (km2)	Ground Deformation (mm/a)		GWD (m)		Change Rate of GWD (m/a)		Compressible Thickness (m)
			Mean	Range	Mean	Range	Mean	Range	Mean	Range
A	I & II	4036	−22.60	−88.30–10.35	24.96	4.19–41.90	−0.13	−1.16–0.66	75	51.40–112.45
	III				46.81	31.58–63.46	−0.28	−1.87–0.60	106	68.62–135.74
	IV				49.91	38.28–68.02	−0.84	−2.48–−0.08	104	78.91–123.96
B	I & II	5746	−23.86	−78.83–4.30	13.43	2.79–36.81	−0.57	−0.96–−0.23	83	54.20–122.68
	III				74.72	52.28–93.48	−1.16	−2.52–0.30	120	67.43–173.67
	IV				79.24	53.77–98.68	−1.30	−2.46–−0.19	95	68.97–129.29
C	I & II	139	−17.45	−46.95 ~−3.78	4.79	4.30–5.32	−0.26	−0.37–−0.22	84	81.77–87.19
	III				69.35	66.32–72.03	−0.38	−1.04–−0.04	135	123.98–147.21
	IV				61.90	59.37–65.70	−0.60	−1.06–−0.27	87	83.14–91.42
D	I & II	675	14.47	7.13–38.30	6.69	4.30–5.31	−0.78	−1.14–−0.55	100	91.01–103.94
	III				49.41	38.99–54.77	−1.84	−2.99–−0.75	138	129.41–160.18
	IV				50.72	46.57–56.61	−1.11	−1.47–−0.89	115	107.68–122.03
E	I & II	832	16.76	0.12–28.27	4.51	2.41–5.73	−0.32	−0.46–0	92	76.21–100.84
	III				62.62	52.25–73.54	−1.22	−2.57–−0.60	84	60.65–104.52
	IV				61.09	53.24–77.13	−1.31	−2.36–−0.40	85	75.94–94.57
F	I & II	716	11.15	4.89–48.36	8.26	6.64–12.78	−0.59	−0.84–−0.39	61	44.01–96.27
	III				65.48	44.30–75.44	−2.38	−3.05–−1.19	69	45.82–149.81
	IV				72.04	65.42–79.34	−2.44	−2.92–−1.70	99	80.39–111.69

), a random dataset comprising 100,000 data samples for each variable ($x_i$) was generated. Subsequently, the predicted deformation values ($Y$) for each sample were calculated using Equation (1). Given a deformation threshold range of −10 mm/a to 10 mm/a, those samples whose predicted values fall within this threshold were then compiled into a dataset.

In the field of geoscience, the mean value or average is frequently employed as a threshold due to its ability to encapsulate the central tendency of data, offering a pivotal reference point for analysis [39,40]. In addition, the 25% and 75% quartiles are commonly employed as threshold ranges for controlling variables in datasets, due to their effectiveness in illustrating the distribution and essential features of the data [41,42]. Thus, in this study, the compiled dataset underwent statistical analysis to determine key metrics: the mean, the 25% quartile, and the 75% quartile of each influencing factor. These values were then established as the thresholds for each influencing factor. Utilizing these thresholds can guide the range of influencing factors, which in turn could potentially aid in controlling ground deformation.

3. Results

3.1. Spatiotemporal Variations of Land Deformation and the Identified Subsidence and Uplift Areas

The mean of annual land deformation rate derived from SBAS-InSAR in the study area is shown in Figure 3. The SBAS-InSAR-derived deformation has been validated using the monthly GNSS data [27,43].

The mean annual land deformation rate during 2017–2022 varies greatly in the whole region (−93–47 mm/a). Generally, a large portion of the region showed subsidence while the uplift was also observed in some of the areas. Here, the deformation rates $\le$−10 mm/a and $\ge$10 mm/a were set as the subsidence and uplift, respectively. This threshold rates agree with the global subsidence classification, with the rate <−10 mm/a considered as the nominal or zero subsidence class [5]. Then we can identify the remarkable subsidence and uplift areas where the gridded boundary values are equal to −10 or 10 mm/a, respectively. The identified subsidence areas are labeled as A, B, and C, and the uplift areas are labeled as D, E, and F in Figure 3. In A, B, and C, the subsidence area is 4036 km², 5746 km², and 139 km², and the areal mean subsidence rate is 22.60 mm/a, 23.86 mm/a, and 17.45 mm/a, respectively (Table 1). Clearly, B has the maximum subsidence area and rate. Moreover, within each of the subsidence areas, the deformation rate changes greatly in space, e.g., −96.20–19.26 mm/a in A, −73.74–21.02 mm/a in B, and −24.30–6.94 mm/a in C. Therefore, it is worth clarifying that a small part of the ground rebound exists within the subsidence area circled using the gridded boundary of −10 mm/a.

In D, E, and F, the uplift area is 675, 832, and 716 km², and the areal mean uplift rate is 14.47 mm/a, 16.76 mm/a, and 11.15 mm/a, respectively. The spatial uplift varies in a range of 2.96–38.21 mm/a in D, 6.64–49.98 mm/a in E, and 2.42–58.47 mm/a in F (Table 1). Clearly, the uplift area and annual change rate are much smaller than those of the subsidence areas.

The cumulative and annual deformation values from 2017 to 2022 at the subsidence areas of A, B, and C, and the uplift areas of D, E, and F are shown in Figure 4. Generally, the subsidence tends to be mitigated, particularly in A and B (Figure 4a,b). In A, the areal mean subsidence rate reduced from 36.35 mm/a in 2017 to 23.27 mm/a in 2022. In B, the areal mean subsidence rate enhanced in the first three years of 2017~2019 and then rapidly mitigated during 2020–2022, resulting in reduction in the subsidence rate in 2022 (13.70 mm/a), compared with 21.91 mm/a in 2017. The annual subsidence rate in C shows a great variation in the period of 2017–2022 (Figure 4c).

In the uplift areas, the cumulative deformation from 2017 to 2022 reached 144 mm in D, 82 mm in F, and 152 mm in E. The annual uplift rate shows a sustained increase in F and a decrease in E during 2017–2022. In D, the uplift rate increased to the maximum of 29.56 mm/a in 2020 and then decreased. Such differences in annual subsidence and uplift rates in these areas may result from different variations in GWDs and deposited properties in space (Table 1).

3.2. Spatiotemporal Variations of the Influencing Factors

3.2.1. Groundwater Table Depth and Its Change Rate

As shown in Figure 5, the three uplift areas (D, E, and F) are located at the areas with a high groundwater level (e.g., D) and the depression cone boundaries in the deep aquifers (III and IV) (e.g., E and F). Therefore, we postulate that the ground rebound could be sensitive to changes in the groundwater table, since the shallow groundwater areas and boundary areas easily produce great variations in the groundwater table and high hydraulic gradients when groundwater extraction reduces or ceases.

The change rate of the annual GWD in space is illustrated in Figure 6. The negative change rate occupied 90.45%, 85.78%, and 97.91% of the study area in the aquifers I & II, III, and IV, respectively, which indicates that the groundwater level has risen widely since 2017. The obvious rise is especially marked for the deep aquifers (III and IV) and in the uplift areas of D, E, and F (Figure 6b,c), with the rising rate of 0.61–3.05 m/a for aquifer III, and 0.41–2.91 m/a for aquifer IV. This indicates that the uplift pattern should be attributed to the head rise due to the restriction of groundwater exploration policy.

The decline in the groundwater level with a maximum rate of 0.60–0.72 m/a is only noticeable in the subsidence area of A for the aquifers I & II (Figure 6a), indicating that the subsidence is mainly due to groundwater declines for A. For the other subsidence areas of B and C, the groundwater level shows an increasing trend except for a small portion of B (2.78% of the B area). This indicates that the subsidence is continuing even though the groundwater level has stopped falling. However, the rising of groundwater level has mitigated subsidence since its annual subsidence rate has reduced (Figure 4).

3.2.2. Thickness of Compressible Sediments

The spatial distribution of thickness for compressible sediments in the three aquifers is shown in Figure 7. For the shallow aquifers of I & II in Figure 7a, the compressible thickness ranges from 43 to 124 m, and is smaller at the southern part of A and covers the largest area in F. In the deep aquifers of III and IV, the small thickness is located at E and F, while the large thickness of over 100 m is located in the areas of D and C, the depression cone in B, and the northern part of A. The proportion of compressible sediments in the whole aquifer in Figure 7d highlights less compressible silt and clay deposits (<50% of the total thickness) in the southern part of A, and the whole area of F and E, and greater compressible sediments (>55% of the total thickness) in the northern part of A, and most of the areas of D, B, and C.

The spatial distribution of compressible sediments indicates the potential for highly elastic (recoverable) areas in E and F, and highly inelastic (unrecoverable) compaction areas in D, B, and C. Thus, rising in the groundwater table should lead to a greater ground rebound in E and F, and less rebound in D, B, and C. However, the high water level in D (<50 m depth) in Figure 5 could increase the pore pressure, thus lowering the effective stress and inelastic compaction. As a result, the ground rebound in D is similar to that in E and F when the groundwater table rises (Figure 6). In contrast, in spite of the rise in the groundwater level, the highly inelastic compaction plus the lowest water table (>70 m depth) in the deep confined aquifers of III and IV leads to significant subsidence in B and C (Figure 5b,c). This indicates that the groundwater level still stays below the threshold for initializing the upward rebound of strata. In A, the rapid decline in the groundwater table plus the large compressible thickness in the north can strengthen the effective inelastic compaction, whereas the rise in the water table and less compressible sediments in the south of A significantly reduce the effective inelastic compaction even though the groundwater table is relatively lower.

3.3. Prediction of Land Subsidence and Uplift

3.3.1. The RF Prediction Accuracy

The training dataset was randomly sampled and thus each training result could be different. In this study, we found that repeating the modeling process ten times using RF can yield stable results (stable R², RMSE, and MAE). The fitted lines between the observed and stably predicted value of deformation in the test set using the least square regression (red line in Figure 8) are close to the 1:1 line (black dotted line in Figure 8). Specifically, R² is 0.95 and 0.93 for the validation dataset in the subsidence and uplift areas, respectively.

A kernel-density estimate is used to represent the correspondence between the observed value and the predicted value (Figure 8). The high frequency (>0.01) is concentrated in the range where the observed and predicted values are nearly equal, which occupies 69.63% and 80.89% of the entire dataset in the subsidence and uplift areas (Figure 8a,b), respectively. In addition, the MAEs are 12.77 mm and 6.03 mm, and the RMSEs are 17.45 mm and 8.65 mm, in the subsidence and uplift areas, respectively. Overall, the RF can well explain the relationship of the gridded land deformation with the driving factors in subsidence and uplift areas.

3.3.2. The Relative Importance of the Influencing Factors

Equation (6) was used to calculate the importance score of each input variable. As a result, the importance scores of features are ranked in Figure 9. In the subsidence areas A, B, and C, the compressible thickness of the confined aquifer III (X8) is highlighted to be the main driving factor for the land subsidence, with the highest score of 21.73% (Figure 9a). The change rate in annual GWD for the confined aquifer IV (X6) ranks second with the importance score of 18.62%. The other factors are the change rate of annual GWD for the confined aquifer III (X4) and the shallow aquifers I & II (X2), as well as the mean annual depth in the shallow aquifer (X1) and the two deep aquifers (X3 and X5), with the scores of 8.16–10.21%. The low-influence variable is the compressible thickness of the shallow aquifers and the deep aquifer IV (X7 and X9), which have the score of 7.75% and 6.94%, respectively.

In the uplift areas of D, E, and F, the highest score (40.04%) is found to be the change rate of annual GWD in the shallow aquifers I & II (X2) (Figure 9b), demonstrating that the rise in the shallow groundwater table contributes most to the ground rebound. Moreover, the change rate of annual GWD in the deep aquifer IV (X6) ranks second with the score of 20.44%. The four other features (X5, X1, X4, and X3) are associated with mean GWD and its change rate, which have scores in the range of 5.19–9.73%.

Among these influencing factors in the subsidence and uplift areas, the annual change rate of GWD in all aquifers contributes a total score of 38.19% in the subsidence area, much smaller than the total score of 65.92% in the uplift area. In particular, the change rate of annual GWD in the deep aquifer IV (X6) ranks second in both subsidence and uplift areas. The importance score of mean annual GWD is almost the same in the three aquifers (8.16–8.63% for X1, X3, and X5 in Figure 9a) in the subsidence area, with a total score of 25.39%. This total score is a little larger than the score of 22.69% in the uplift area. In contrast, the importance of compressible thickness of all aquifers contributes a total score of 36.42% in the subsidence area, much higher than the value of 11.38% in the uplift area. The thickness in the confined aquifer III (X8) is emphasized by comparing the scores with other aquifers. These results suggest that the change in the groundwater level plays a significant role in land deformation, particularly for ground uplift when the groundwater level rises in the shallow aquifers (e.g., X2). The changes in the groundwater level have almost the same effect as the compressible thicknesses on land subsidence (e.g., X6 and X8).

3.3.3. The Appropriate Values of Influencing Factors for Mitigating Land Deformation

It is well known that land subsidence has caused multiple adverse impacts, while fast rebounding could also result in geohazards. For example, a fast recovery of the piezometric head caused infrastructure damages by buoyant forces acting on the building foundations and groundwater seeped into the basement floor of buildings and tunnels in Tokyo [44]. Therefore, the subsidence and uplift should be mitigated via controlling their influencing variables, such as the groundwater level and its change rate. In our study region, the deformation rate within −10–10 mm/a is widespread and has yet to cause severe environmental disasters. Thus, from the RF generations of deformation rates in relation to possible changes in the driving factors, we selected the gridded values of the combinations of influence variables that meet the deformation rate within −10–10 mm/a (Figure 10). Figure 10 and

Table 2. The measured and controlled mean GWDs and change rates in different aquifers and areas (note: M refers to the measured value, and C refers to the controlled value).

Area	Mean GWD (MGD, m)						Change Rate of MGD (m/a)
	I and II		III		IV		I and II		III		IV
	M	C	M	C	M	C	M	C	M	C	M	C
A	24.96	15.22 (7.04–19.7)	46.81	53.59 (41.72–64.27)	49.91	58.94 (49.22–67.75)	−0.13	−0.69 (0.59–0.77)	−0.28	−1.50 (1.34–1.65)	−0.84	−1.53 (1.05–2.06)
B	13.43		74.72		79.24		−0.57		−1.16		−1.30
C	4.79		69.35		61.90		−0.26		−0.38		−0.60
D	6.69	4.80 (3.48–5.69)	49.41	54.94 (42.90–66.91)	50.72	70.79 (67.17–75.31)	−0.78	−0.15 (0.05 ~0.23)	−1.84	−1.10 (−0.64–−1.58)	−1.11	−1.45 (−0.54–−2.27)
E	4.51		62.62		61.09		−0.32		−1.22		−1.31
F	8.26		65.48		72.04		−0.59		−2.38		−2.44

exhibit that, in the subsided areas of the shallow aquifers I & II, the areal mean GWD (MGD) should be controlled to 15.22 m (ranging within 7.04 m–19.73 m for the 25% and 75% quartile) and its elevated rate is controlled to be 0.69 m/a (range of 0.59–0.77 m/a), respectively. In aquifer III, MGD and its elevated rate should be controlled to be 53.59 m (range of 41.72–64.27 m) and 1.50 m/a (range 1.34–1.65 m/a), respectively. In aquifer IV, MGD and its elevated rate should be controlled to be 58.94 m (range of 49.22–67.75 m) and 1.53 m/a (range of 1.05–2.06 m/a), respectively. In particular, compared with the measured MGD and elevated rate (Table 2), the controlled MGD should be smaller than 19.7, 64.27, and 67.75 m (the 75% quartile in Figure 10) in the depression cone centers for the aquifers of I & II, III, and IV, respectively. Meanwhile, the rising rate of the groundwater level should be increased to 2–5 times of the current rates in the depression cones of A and C. Moreover, as the compressible thickness of the confined aquifer III (X8) is the most important driving factor, a compressible thickness of less than 84.91 m (65.09 m–104.36 m) should be retained in the aquifer subsidence (Figure 10).

The uplift areas of E, F, and D are located at the shallow aquifers of I & II and around depression cones of the deep aquifers of III and IV (Figure 5). In the shallow aquifers of I and II, MGD is suggested to be 4.8 m (3.48 m–5.69 m) and maintains a slowly rising rate of 0.15 m/a (0.05–0.23 m/a). This controlled rate is about 1/2–1/5 of the measured rates in D, E, and F (Table 2), indicating a slowly rising rate can mitigate the uplift. In the deep aquifers, reducing the rise rate of the groundwater table, particularly for F (e.g., −1.0 m/a), can prevent the significant uplift.

4. Discussion

4.1. The Effect of Groundwater Level Rising on Land Subsidence

The rising of the groundwater level has been widely distributed in our study region, especially in the deep aquifers in recent years (Figure 6). However, large areal ground subsidence still occurred, such as in B and C (Figure 3). Holzer [45] suggested that significant subsidence in the field might not be observed until hydraulic heads declined >30 m in many aquifer systems, or until the natural pre-consolidation stress is reached. Our analysis shows that the three remarkable subsidence areas (A, B, and C) are located at groundwater depression cones where groundwater table decreases by >40 m in the shallow aquifers for A, and >90 m in the deep aquifers for B and C (Figure 5). Moreover, the highest subsidence rate of over 90 mm/a in the depression cone is much larger than the historical highest value of over 50 mm/a found in Cangzhou during 1971–2015 [46], and 60 mm/a during 1990–2010s in the northwest of Cangzhou [25,47]. Moreover, the regions with high compressible thickness (especially higher than 84.91m in the third aquifer) are more likely to experience land subsidence, which should be noticed. This suggests that the rise in the groundwater table has not raised the effective stress to be higher than the natural pre-consolidation stress after the prolonged water table decline in the study region. Compared to sand with hard quartz grains, a large proportion of clay in our study areas is prone to compaction when increased load upon it or reduced support within it can cause its property of being easily deformed [48].

Nevertheless, the subsidence rate tends to be reduced in some areas as the groundwater table rises (Figure 4). For example, the areal mean subsidence rate in A reduced from 36.35 mm/a to 23.27 mm/a during 2017–2022 because of ceasing or reducing groundwater abstraction in recent years. Here, in terms of the gridded mean values of subsidence accounted in an interval (e.g., 0.1 m/a in Figure 11a) of the change range in the subsidence areas, we illustrate the relationship of the subsidence rate with each of the influencing factors in A, B, and C (Figure 11). Figure 11 clearly shows that annual land subsidence is generally proportion to the compressible thickness in deep aquifers (Figure 11a), and the groundwater table depth in the aquifers (Figure 11b,c,e). So, the increase in the groundwater level and consequent reduction in the compressibility of the clay could substantially mitigate subsidence.

Figure 11 also exhibits that the multidimensional data have nonlinear correlation structures between the predictor (land subsidence) and variables (driving factors). As shown in Figure 11a, the subsidence sharply increases with the compressible thickness of confined aquifer III when the thickness is smaller than 120 m and then fluctuates around a steady value.

4.2. The Effect of Groundwater Table Rising on Land Uplift

An appreciable portion of the study region shows ground uplift. The largest areal mean uplift rate is 16.76 mm/a in E and could reach over 50 mm/a at some specific sites of the area. With the rising groundwater level, uplift tends to happen because of the related effective stress unloading. This uplift rate was larger than 12.6 mm/a when the over-exploitation of groundwater was initially banned and the hydraulic head rose quickly to 37.6 m in Cangzhou during the early 2000s (2000–2013) [28]. The largest uplift rate of over 50 mm/a in the depression cone center is also larger than the uplift center rates of up to 25 mm/a during 2017–2020 in the Taiyuan basin of North China [49]. In the aquifer characterized by multiphase materials, the coupling between pore fluid and the soil skeleton can be considered elastic [16]. The higher uplift rate with the rise in the groundwater level in recent years suggests that the consecutive elevation in the groundwater level in the stratified aquifer systems could stimulate an increase in elastic porous storage and aquitard re-pressurization, leading to the poroelastic rebound mostly in the shallow aquifers. To control subsidence, MGD less than 20 and 70 m is suitable in the shallow and deep aquifers, respectively, and the increasing rate of MGD should be higher than 0.69 m/a and 1.50 m/a in the shallow and deep aquifers, respectively.

The relationship of land uplift rate with each of the influencing factors in the uplifted areas (D, E, and F) is illustrated in Figure 12. The uplift rate is proportional to the negative change rate of GWD, particularly in aquifers of I & II and IV, suggesting that the rise in the groundwater level in aquifers can significantly make the ground rebound. The nonlinear correlation structures between the dependent and independent variables also exist, as shown in Figure 12. In Figure 12a, when the shallow groundwater level rises at a rate larger than 0.5 m/a, the rebound rate reaches its peak value and does not respond to the further change in the groundwater level. In Figure 12b, the effective ground rebound can be found only when the rise rate of the groundwater table in deep aquifers is larger than 2.0–2.4 m/a. This indicates that a smaller rise in the groundwater level in aquifer IV may not significantly elevate pore pressure to a critical value (e.g., 60 m) that leads to a ground rebound (see Figure 12b).

Comparatively, the uplift areas and rates are smaller than those of subsidence since clay subsidence is largely a one-way process as its inelastic property cannot be recovered [49,50]. Thus, land subsidence can permanently reduce available groundwater, especially in zones with large compressible thickness. The estimated permanent loss of groundwater storage capacity can reach 9.41 billion m³ per meter decline of the confined water level in Cangzhou [18]. The inelastic compaction accounted for 87.4–87.9% of total groundwater storage depletion in Cangzhou [24].

5. Conclusions

In this study, SBAS-InSAR data from January 2017 to December 2022 were used to delineate the land deformation in the multi-layered porous aquifer system in the Cangzhou and Hengshui region of the NCP. We found that even though the groundwater table has risen widely in the study area, land subsidence and uplift, with a rate of >10 mm/a, co-exist. In the identified three remarkable subsidence areas, the areal mean subsidence rate ranges from 17.45–23.86 mm/a, mostly located in the depression cones where the depth to the GWD has lowered to over 40 m and over 90 m in the shallow and deep aquifers, respectively. Even though the subsidence rate significantly reduces with rising of groundwater level, the inelastic compaction continues, as shown in the accumulated subsidence, indicating that effective stress is still higher than the natural pre-consolidation stress. Meanwhile, three noticeable land uplift areas are identified. The uplift rate ranges from 11.15 to 16.76 mm/a, mostly located in the declining areas of the groundwater level in the shallow aquifers and the depression cone boundaries in the deep aquifers where the groundwater table rose markedly during the study period.

The RF method was proven to be effective for identifying complex relationships with high dimensionality between land deformation and the influence variables, such as mean GWD and change rate, as well as compressible thickness in each aquifer. R² is 0.93 for the validation dataset in the subsidence areas and 0.95 in the uplift areas. The importance scores based on the trained RF successfully delineate the importance of the influencing factors in the subsidence and uplift areas. Since the subsidence occurs in the depression cones, the thick compressible layer of the confined aquifer is the main driving factor, whereas the rise in the GWD plays an important role in reducing the inelastic compaction. Comparatively, the fast rises in the groundwater table in the shallow aquifers (−0.34 m/a) and the deep aquifers (−1.14 m/a) are the main driving forces for making ground elastically rebound.

Mitigation of the fast land deformation, including subsidence and uplift, could prevent occurrence of environmental hazards. Based on the RF predicted possible deformation rates resulting from changes in the driving factors, we selected appropriate combinations of driving factors that could hamper ground deformation rates. To control the subsidence, the appropriate MGD is less than 20 and 70 m in the shallow and deep aquifers, respectively, and the rising rate of MGD increases to 2–5 times of the current rates in the depression cones. The exerting pressures in the confined aquifers could avoid compressibility of the clay for the compressible thickness <84.91 m. To mitigate the uplift, the appropriate rising rate should be reduced to 1/2–1/5 of the current rates in the shallow aquifers.

The land subsidence can permanently reduce available groundwater, especially in zones with a large compressible thickness. Uneven land deformations, such as subsidence at the depression cone center and uplift at its boundary, may exacerbate geohazards. Therefore, it is vital to implement appropriate governance of groundwater recovery for hazard prevention. Furthermore, in this study, the data of groundwater level are annual; the research could be improved by using monthly data for the future study of seasonal features or the time lag between them and land deformation under the groundwater exploitation restriction policy.