The study area is in the upper and middle part of the Chaobai River’s alluvial fan, which was formed in the Pleistocene by the alluvial-diluvial action of the Chaobai and Wenyu Rivers (Lei et al., 2022a). The geological environment in this area is highly complex due to the heterogeneous distribution of alluvial deposits (Wei, 2008). The thickness of the Quaternary deposits varies between 50 m in the northeastern region and extends up to 1,100 m in the southwestern region, change from an coarse sand in the piedmont areas to fine clay in the floodplain zones (Zhu et al., 2020).

Two buried geological faults are located: the northern section of the Gaoliying (GLY) fault (Figure 1A) and the northern section of the Shunyi (NSY) fault (Figure 1A). The GLY fault is a tension normal fault with a length of 40 km, a strike of 70°–80°, and a SE trend (Wei, 2008). The GLY fault has experienced multiple periods of activity since the Cretaceous Period and has been very active since the Cenozoic in controlling the Palaeocene, Neoproterozoic, and Quaternary deposits (Wei, 2008). The average activity rates of the Early Pleistocene, Middle Pleistocene, Late Pleistocene and Holocene faults were 0.07 mm/a, 0.04 mm/a, 0.23 mm/a and 0.10 mm/a, respectively (Zhang et al., 2016).

The NSY fault is a tension normal fault with a length of 40 km, a strike of 40° and a SE trend (Wei, 2008). According to the measurements and geochronological results, the average activity rates of the Early Pleistocene, Middle Pleistocene, Late Pleistocene and Holocene faults were 0.23 mm/a, 0.03 mm/a, 0.29 mm/a and 0.51 mm/a, respectively (Qi et al., 2020).

The aquifers in the study area can be divided into four groups (Figure 1B). Phreatic aquifers (PAs), the first confined aquifers (FCAs), the second confined aquifer (SCA), and the third confined aquifer (TCA). Phreatic aquifers (PAs) originated in the Holocene and late Pleistocene and are widely distributed in the study area with thicknesses varying from 40 m to 50 m. Medium-coarse sands with gravel were deposited along the PA layers. The PA is supplied by atmospheric precipitation, infiltration from agricultural irrigation and river water, and its own abundant groundwater resources. The first confined aquifers (FCAs) were formed in the mid-Pleistocene, with depths ranging from 80 to 120 m. The second confined aquifer (SCA) is mainly distributed in the middle of the alluvial fan. Its layers were formed in the lower Pleistocene. The bottom depth of this group ranges from 150 to 180 m. The third confined aquifer (TCA) is distributed south of Mapo village. It was formed in the Lower Pleistocene, and the bedrock comprises the bottom of this group.

The MAR commenced in 2015, with recharge operations conducted annually from April to June. By the end of 2020, a total of 501, 508, 900 cubic meters of groundwater had been recharged. The South-to-North Water Diversion water, after being released through the Lishishan Diversion Gate on the Jingmi Diversion Canal, flows through the Xiaozhong River, Mangniu River, and Huaihe River before entering the upstream groundwater recharge site at the Chaobai River Niulanshan Rubber Dam (Lei et al., 2022a).

Since its launch in 2007, the Radarsat-2 satellite has been in operation. With a revisit period of 24 days, it employs a C-band Synthetic Aperture Radar (SAR) featuring a 56 mm wavelength. The satellite boasts a spatial resolution of 30 m, making the data it collects ideal for monitoring regional subsidence (Ng et al., 2017; Samsonov et al., 2017).

In this study, 70 scenes of the Radarsat-2 dataset in W1 mode between October 2014 and October 2021 were collected to acquire the regional subsidence of the study area. Persistent scatterer InSAR (PS-InSAR) technology was adopted to process the SAR dataset and extract subsidence information. The basic principle of the PS-InSAR technique is to use several SAR images of the same area to find permanent scatterers (such as buildings, bridges, roads, etc.,) that are not subject to temporal and spatial baseline deconvolution and atmospheric effects by statistically analysing the magnitude information of the images and then to separate the errors from the differential interference phases one by one using the network formed by connecting these PS points and the spatiotemporal characteristics of each phase component (Ferretti et al., 2001). Using the network of these PS points and the spatial and temporal characteristics of each phase component, the errors are separated from the differential interference phase one by one, and the surface deformation phase of each PS point is finally obtained (Ferretti et al., 2001). It can overcome the factors of spatial and temporal incoherence and atmospheric delay in conventional differential interferometry and thus obtain continuous and reliable surface deformation information (Ferretti et al., 2001). The accuracy of the annual average deformation rate obtained by the PS-InSAR method can reach the millimetre level (Lei et al., 2022b). The data were processed via the following steps:

The master image is determined based on information such as the spatial and temporal baseline, the Doppler shift and the number of datasets. The remote sensing images are then cropped according to the extent of the study area to determine the processing area of the remote sensing images. The coordinate mapping relationship between the master image and all SAR imagery in the distance and azimuth directions is calculated, and the mapping relationship is used to perform coordinate mapping on all SAR imagery. Then, the image is resampled to a new image to align all the SAR images to the master image. The 30 m resolution SRTM data are resampled into the SAR image coordinate system, and an amplitude image is generated from the reference SRTM data. The aligned images were then subjected to differential interferometry to obtain N differential interferograms, with each image element of each differential interferogram containing five components. To avoid calculation errors caused by low coherence points, the high coherence points in the differential interferograms were selected for calculation. In this paper, the amplitude deviation index is used for screening, and the time series amplitude extreme deviation values are obtained by calculating the mean and standard deviation of the PS points in each period and the ratio of the two. The threshold value of the amplitude deviation index is set to filter the PS points, and then the Delaunay triangulation network is used to construct the point network of the PS point analysis. Once the PS points are selected, strong reflection points with zero phase components of deformation and elevation, high coherence and stable reflection characteristics are selected as the reference points. Then, the phase distributions of all other pixels are calculated based on the current reference points, and the phase values are improved based on the above phase distribution. The surface deformation of the study area and the physical and statistical properties of the components of the interferometric phase are used to establish an entangled phase function model with the elevation error, linear deformation, thermal deformation and other factor parameters. The atmospheric phase is analysed according to the pixel network using the phase difference relationship between adjacent permanent scattering targets, and the function model is solved based on the inverse residuals algorithm. The phase of the elevation error, linear deformation and other contributing factors of each pixel is solved, the phase of nonlinear deformation and incoherent noise in the residuals is decomposed, and finally, the entangled phase values of atmospheric delay components such as the elevation difference and linear deformation rate are accurately calculated. The above steps were completed using the SARPROZ tool.

The J1 land subsidence monitoring station is equipped with extensometers of varying depths. Data from these extensometers, spanning from 2004 to 2021, have been selected for analysis to discern the characteristics of soil deformation and to explore the correlation between the deformation of stratified soil layers and groundwater levels. The extensometer data at varying depths can be analyzed by deducting the cumulative deformation of the underlying layer from the total deformation. The monitoring layers of the extensometers and the corresponding groundwater wells at the monitoring station are detailed in Table 1. The geographical location of the monitoring station is depicted in Figure 1.

Data collected from monitoring wells in groundwater level directly reflect the hydraulic parameters of the aquifer. In this study, GWL was collected from 53 monitoring wells. The data is gathered on a monthly basis by the Beijing Institute of Geo-Environmental Monitoring from 2014 to 2021. In addition, time series analysis was conducted using data from 6 observation wells at both sides of the NSY fault. W1, W2, S1, S2, S3 and S4 wells. The depths of the wells are 111 m, 102 m, 60 m, 50 m, 80 m, and 102 m, respectively. In this study, the groundwater level refers to the elevation of groundwater, and the calculation method is to subtract the groundwater depth from the surface elevation.

Fifteen levelling benchmarks evenly distributed across the study regions were measured to verify the reliability of the PS-InSAR processing results (Figure 2A). Owing to the absence of permanent scatterers (PS) at the reference points, the average displacement of PS points within a 100-m buffer zone surrounding these leveling points was employed as the land deformation value in the reference point during the calibration process.

As shown in Figure 2B, a cross-comparison levelling result in 2014 and 2021 and the PS-InSAR data was obtained via linear regression analysis. The assessment results reveal that the PS-InSAR data matched the levelling results, and the error ranged from 1 to 7.5 mm. The R-squared value was 0.976. The comparison validated the PS-InSAR data, and the data could be used for further analysis of the spatial pattern of land deformation.

As depicted in Figure 3, a comparison of the velocity maps from 2014 to 2021 reveals that both land subsidence relief have been observed within the study area. The deformation values higher than 0 indicate uplift, while the values lower than 0 indicate land subsidence.

Significant land subsidence has been detected in the southern region of the study area, with a maximum deformation rate of −83.8 mm/yr. Conversely, in 2021, widespread uplift was documented in the Shunyi district, with rates varying between +10 and +20 mm/yr (Figure 3). The PS (Persistent Scatterer) data within the land subsidence funnel indicate subsidence rates exceeding −50 mm/yr in 2014. By 2021, the region experienced land uplift at a rate as high as +5 cm/yr. The maximum land rebound was observed in the region around 2021, amounting to ten percent of the land subsidence recorded in 2014.

Interesting, the visible uplift zone, the rebound rate above 20 mm/yr, is dominated spreading in a fan area striking approximately NE‒SW. By jointly preexisting faults and land deformation driven by the InSAR, the uplift zone is mainly developed in the depression bounded by the GLY fault in the northwest and the NSY fault in the southeast.

For further reveal the effect of pre-existing fault on the land deformation, We drew a profile a-a (marked as the brown dashed line in Figure 3B) perpendicular to the pre-existing fault. We plotted the trend of vertical displacement (black line) and the gradient of subsidence (brown line) based on the points (600 points total) within the profile. As shown in Figure 3C, the deformation rate reveals an decreasing trend at the hangingwall of NSY and footwall of GLY. The maximum value of the vertical displacementis 38.89 mm. A visible land rebound area were generated at the foot wall of NSY, with a maximum rebound of 18.54 mm. The deformation rate exhibits a distinct uneven changes crossing the GLY Fault and NSY fault, with a maximum uneven deformation of 18.7 mm.

Previous studies have shown that land uplift could result from the expansion of the granular skeleton due to the increase in pore pressure associated with recovery of the groundwater level (Waltham, 2002; Wang et al., 2017). This phenomenon is referred to as elastic or poroelastic rebound (Chen et al., 2007). Therefore, we consider the change in groundwater levels on both sides of the preexisting faults to likely be triggering factors that influence the spatial extent of land rebound.

53 automatically monitored wells were used to investigate the groundwater level changes from 2014 to 2021 at both sides of the preexisting faults. As shown in Figure. During the MAR operation, Groundwater levels exhibited a relative increase at 90% of the monitoring well sites. Groundwater levels in the footwall of the NSY fault have been observed to rise, with a peak increase of nearly 50 m. In contrast, the groundwater levels at the hanging wall of the NSY fault exhibited a relatively minor upward trend. As the distance from the fault increased, no discernible changes in groundwater levels was observed.

To further examine groundwater level changes at both sides of the preexisting fault, The monitoring wells, located at at both sides of the fault, was selected. The location of the wells are shown in Figures 4, 6 respectively. The groundwater level exhibited significant difference on either side of the NSY fault. For instance, monitoring well 3 showed an visible increase in response to the MAR with 10.71 increase after the MAR operation in 2018. While W1 record seasonal variations associated with changes in precipitation throughout the observation period Figure 5A.

A similar finding was made by Lapperre et al. (2022) for the Peel Boundary Fault. They concluded that this effect may caused by low normal permeability of the preexisting fault zone. Previous studies have shown that the most direct cause of the low permeability around the fault zone is because of fault throwing, which can juxtapose low-permeability confining units against well-permeability aquifers (Lapperre et al., 2019; Zhao et al., 2021).

The geological section of the NSY fault at a depth of 0–100 m was revealed by drilling. As shown in Fig, notable stratigraphic dislocation begins at a depth of 35 m, and the vertical throw of the strata gradually increases, displaying a juxtapose of low-permeability confining units against well-permeability aquifers (Figure 6). The geophysical prospecting profiles also show that the stratigraphic dislocation propagates downward into the sediments and may extend to depths of more than 100 m (Qi et al., 2020).

The earth fissure trenches reveal the shallow characteristics of pre-existing faults. There is a clear vertical offset of strata within the trench’s range, resulting in the development of a clay smear along the fault zone that cuts through the clay beds (Figure 7).

Based on the results of the trenching and geological section analysis, significant stratigraphic displacement and distinct fault surfaces were observed. The results indicates that the NSY Fault displays pronounced brittle fault characteristics. As shown in Fig, under the effect of the MAR, the water level rise was occurred in the foot wall of NSY fault Figure 5B. The juxtapose of low-permeability confining units against well-permeability aquifers caused by the throwing lead to a gradual decrease in the hydraulic conductivity of the fault zone. An obvious reduction in the permeability of half the magnitude in the hanging wall of the NSY fault zone compared to the aquifer in the footwall of the NSY fault zone was found (Table 2).

As shown in Figure 8B The presence of a low-permeability layer diminishes the hydraulic connection between the hanging wall and footwall of the NSY fault, the groundwater in the deep aquifers located in the footwall of the NSY fault have to circumvent the aquiclude by moving upward and downward through aquifers situated above and below it,resulting in reduced water transmission from the MAR to the hanging wall of the NSY fault, as compared to the natural state Figure 8A. Consequently, this results in a hysteresis effect or a lack of groundwater level rise in response to the MAR.

According to the principle of effectiveness, MAR leads to an obvious rebound of the aquifers located at the footwall of the NSY fault. On the other hand, due to the sluggish recharge of water from the footwall to the hangingwall aquifers, there was few rebound in contrast to the deformation observed in the footwall of the NSY fault.

Ground surface deformation is a result of the total deformation of deep soil layers (Liu et al., 2019; Zhang et al., 2015). It can be inferred that the spatial variation may be caused by the complex loading and unloading in deep soil deformation caused by the groundwater level rise. To better understand the relationship between the groundwater level and the deep soil deformation, the deformation of the soil layer measured at extensometer station, located at the foot wall of the NSY fault, is plotted in Figure 9.

As shown in Figure 9, the land rebound responds differently to a rise in groundwater level. Land rebound can be broadly grouped into two deformation modes.

Firstly, as the groundwater level rises, the deformation of the aquifer is closely related to the change in the groundwater level and exhibited the soil rebound. Secondly, the deformation of the aquifer is no closely related to the change in the groundwater level. The relieve of soil compaction was found during the monitoring period following the rise in water levels. Based on the boreholes data, the former mainly occurred in the soil with higher proportion of sand layer, while the latter mainly occurred in the soil with higher proportion of the clay layer.

We made comparisons of the lithology at both sides of the NSY fault to assess the effect of the lithology on the land rebound. As shown in Figure 10, The fairly high proportion of sand layer, suggested that this area may occurred significant elastic or poroelastic rebound in soil during the groundwater level rise. In contrast, lower proportion of sand layer was identified in hanging wall of NSY fault. The information indicates that this area can not experience the rebound immediately when the groundwater level rise.