Factors driving surface deformations in plain area of eastern Zhengzhou City, China

Zi-jun Zhuo, Dun-yu Lv, Shu-ran Meng, Jian-yu Zhang, Song-bo Liu, Cui-ling Wang

1 Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences, Shijiazhuang 050061, China.2 China University of Geosciences (Beijing), Beijing 100083, China.3 Key Laboratory of Quaternary Chronology and Hydro-Environmental Evolution, China Geological Survey (CGS) & Hebei Province, Shijiazhuang 050061, China.

Abstract: With the rapid socio-economic development and urban expansion, land subsidence has emerged as a major environmental issue, impeding the high-quality development of the plain area in eastern Zhengzhou City, Henan Province, China.However, effective prevention and control of land subsidence in this region have been challenging due to the lack of comprehensive surface deformations monitoring and the quantitative analysis of the factors driving these deformations.In order to accurately identify the dominant factor driving surface deformations in the study area, this study utilized the Persistent Scattered Interferometric Synthetic Aperture Radar (PS-InSAR) technique to acquire the spatio-temporal distribution of surface deformations from January 2018 to March 2020.The acquired data was verified using leveling data.Subsequently, GIS spatial analysis was employed to investigate the responses of surface deformations to the driving factors.The findings are as follows: Finally, the geographical detector model was utilized to quantify the contributions of the driving factors and reveal the mechanisms of their interactions.The findings are as follows: (1) Surface deformations in the study area are dominated by land subsidence, concentrated mainly in Zhongmu County, with a deformation rate of -12.5--37.1 mm/a.In contrast, areas experiencing surface uplift are primarily located downtown, with deformation rates ranging from 0 mm to 8 mm; (2)Groundwater level, lithology, and urban construction exhibit strong spatial correlations with cumulative deformation amplitude; (3) Groundwater level of the second aquifer group is the primary driver of spatially stratified heterogeneity in surface deformations, with a contributive degree of 0.5328.The contributive degrees of driving factors are significantly enhanced through interactions.Groundwater level and the cohesive soil thickness in the second aquifer group show the strongest interactions in the study area.Their total contributive degree increases to 0.5722 after interactions, establishing them as the primary factors influencing surface deformation patterns in the study area.The results of this study can provide a theoretical basis and scientific support for precise prevention and control measures against land subsidence in the study area,as well as contributing to research on the underlying mechanisms.

Keywords: PS-InSAR; GIS spatial analysis; Geographical detector model; Degree of contribution of a driving factor; Spatially stratified heterogeneity

Introduction

Surface deformations refer to regional geological phenomena characterized by the gradual uplifting or subsidence of ground elevation, surface inclination, and surface dislocation caused by uneven stress distribution on the surface or underlying soils, influenced by both natural and anthropogenic factors (Xu et al.2009; Chen et al.2023).Excessive groundwater exploitation and large-scale urban construction, driven by urban development have contributed to the global environmental geological issue of land subsidence (Tomás and Herrera, 2010; Yu et al.2020; Yang et al.2022).Regions experiencing rapid economic development in China, such as the Beijing-Tianjin-Hebei,Yangtze-River-Delta, and the Central Plains urban agglomerations, are particularly susceptible to significant urban land subsidence, leading to geological hazards such as surface cracking and collapse, posing severe threat to the safety of urban residents and their properties (Yang et al.2022;Dong et al.2018).Therefore, it is crucial to employ suitable monitoring techniques to obtain information on urban surface deformations, comprehend their spatio-temporal evolution, and systematically investigate their driving factors (Li et al.2018; Liu et al.2022).

Compared with conventional monitoring methods such as leveling measurements and borehole extensometers, the InSAR technique has several advantages.It enables continuous monitoring, irrespective of the time of day or night, and provides high spatial resolution and extensive coverage (Li et al.2017; Liu et al.2018).The PS-InSAR technique, introduced by Ferretti et al.(2000; 2001),can be used to obtain phase information of PS pixels from multiple scenes of time-series SAR images and then perform accurate phase modeling and deformation calculations on the entire scenes of images.This technique overcomes the influence of spatio-temporal decorrelation and atmospheric delay.Concequently, the PS-InSAR technique allows for the timely and precise acquisition of regional time-series deformation information,significantly enhancing the monitoring accuracy for surface deformations.Many domestic and international researchers have utilized PS-InSAR to investigate surface deformation monitoring on various temporal and spatial scales in Chinese cities such as Beijing, Tianjin, and Shanghai,which are affected by severe land subsidence(Chen et al.2016; Liu et al.2016; Dong et al.2018).However, there are limited monitoring results of surface deformations in areas undergoing rapid urbanization, such as the megacity Zhengzhou in the Central Plains Urban Agglomeration, impeding the precise prevention and control measures for regional land subsidence and sustainable urban development.

When examining the factors driving surface deformations, many researchers, both domestically and internationally, have qualitatively analyzed the impact of human activities (e.g.groundwater exploitation and building loads) and natural factors (e.g.hydrogeological conditions) on surface deformations through spatial analysis and mathematical statistics (Ciampalini et al.2016;Zhou et al.2018; Chen et al.2020).Additionally,some scholars argue that the occurrence and development of surface deformations are influenced by multiple factors.For instance, Shi et al.(2018)identified the combined effects of groundwater exploitation and building loads on land subsidence in Shanghai through grey correlation analysis.Cheng et al.(2021) analyzed the spatial heterogeneity of land subsidence in the Chaobai River alluvial fan in Beijing using PS-InSAR and Moran's I.Furthermore, they determined the dominant factor in land subsidence in the middle and lower reaches of the Chaobai River using a geographical detector model.Zhang et al.(2022) investigated the variations in the spatio-temporal patterns of land subsidence during different periods in Wuhan based on a geographical detector model, as well as influencing factors such as hydrogeological conditions, ground load, and underground space development.

The plain area in eastern Zhengzhou covers the main urban area of Zhengzhou and Zhongmu County which has the highest urbanization rate in Zhengzhou.Due to the requirements of urban economic development, the study area has experienced long-term excessive groundwater exploitation, resulting in severe land subsidence (Guo et al.2007; Guan et al.2019; Pan et al.2020).Previous monitoring efforts mainly focused on surface deformations in the main urban area, leaving a limited understanding of the evolution of surface deformations in Zhongmu County.In addition, the analysis of factors driving surface deformations in the study area has been mostly qualitative.As a result, their contributive degrees have not been accurately assessed, impeding precise prevention and control measures for regional land subsidence and hindering sustainable urban development.

Considering the limited research on surface deformations in the study area and the environmental challenges restricting regional sustainable development, this study initially obtained the evolutionary characteristics of surface deformations in the study area from January 2018 to March 2020, which was accomplished by employing the PS-InSAR technique, utilizing 60 scenes of Sentinel-1 images that cover the study area.Subsequently, the study conducted GIS spatial analysis to explore the responses of surface deformations to groundwater exploitation, lithology, and urban construction.Based on this analysis, quantifiable driving factors were identified, and their contributions to surface deformations were quantitatively assessed using the factor detector within the geographical detector model, from which the dominant driving factor was identified.Finally, this study investigated the effects of the interactions among various driving factors on the spatially stratified heterogeneity of surface deformations in the study area using the factor interaction detector.This study aims to provide a theoretical foundation and scientific support for precise prevention and control measures for land subsidence, as well as urban construction planning in the study area.

1 Study area

The plain area in eastern Zhengzhou is located in the southwest of the North China Plain, with geographical coordinates ranging from 113°27′E to 114°12′E and 36°26′N to 35°00′N.It consists of eight administrative districts, including those within the urban area of Zhengzhou (e.g.Huiji District, Jinshui District, Zhongyuan District,Erqi District, Guancheng Hui District, Hightech Industrial Development Zone, and Zhengdong New District) as well as Zhongmu County.The area is divided into three functional planning areas, namely the main urban area, the eastern new urban area, and Airport Economy Zone, covering a total area of 2,195 km2(Fig.1).In recent years,guided by the strategy of accelerating the construction of national central cities, the eastern new urban area has played a critical role in promoting regional economic development and serving as a supplement to the functions of the main urban area.

The study area is dominated by two types of landforms: The alluvial-proluvial inclined plain and the Yellow River alluvial plain.Extensive Quaternary deposits have formed in the region under the influence of neotectonic movements and multiple alluvial processes and flooding events by the Yellow River.These deposits exhibit a gradual increase in thickness and a decrease in grain size from southwest to northeast (Guan et al.2019; Pan et al.2020).The aquifer system in the study area can be divided vertically into two major aquifer groups.The first aquifer group, comprising shallow phreatic water to micro-confined water,consists of Holocene strata encompassing Yellow River alluvium and diluvium, as well as the upper section of Middle Pleistocene strata.This aquifer group is dominated by medium- and finegrained sands.Its footwall is buried at a depth of 45-55 m, under which lie aquitards consisting of silty clay or silty soils, with a thickness of 25-45 m.The second aquifer group, containing moderately deep confined water, is a rock formation buried at depths ranging from 60 m to 360 m.This aquifer group serves as the main source of domestic and industrial water supply in the study area.The second aquifer group can be further divided into two sublayers according to their burial conditions and the physical and mechanical properties of the soils.The first sublayer (II1) is composed of Middle-Lower Pleistocene alluvium,diluvium, and lacustrine deposits.It contains 1-7 aquifers comprising medium- to fine-grained,medium- to fine-grained, and fine-grained sands.The burial depth of the footwall increases from 100 m to 220 m from west to east, with its lower boundary marked by the interface between Quaternary and Neogene strata.The second sublayer (II2)consists of 7-11 beds of medium-grained to coarsegrained, medium-grained to fine-grained, and finegrained sands in the upper member of the Neogene Minghuazhen Formation.The burial depth of its footwall ranges from 167 m to 360 m, gradually increasing from the southwest to the northeast(Fig.2).

Fig.1 Geographic location of the study area and the scope of Sentinel-1A images

2 Data and methods

2.1 Data sources

The data used in this study mainly included Sentinel-1A satellite images, geographic information data, and geological environment data.The Synthetic Aperture Radar (SAR) images were obtained from 60 ascending scenes of images (30 scenes covering the southern and northern parts each) from January 2018 to March 2020, which were acquired in the Interferometric Wide (IW)swath mode from the ESA Sentinel-1A satellite using VV polarization (Fig.1).To remove topographic phase information, elevation data from the SRTM DEM with a resolution of 30 m were employed.This data was sourced from the Resource and Environment Science and Data Center, Chinese Academy of Sciences.Land use data with a resolution of 30 m, were obtained from the long-time-series (1985-2019) annual China land cover dataset (CLCD), and developed by Yang and Huang (2021) from the GeoData Institute.Information pertaining to urban buildings and road networks was derived from the Open Street Map.In addition, the thickness of cohesive soils, the leveling measurements of land subsidence, and the monitoring data of groundwater level were provided by the Institute of Hydrogeology and Environmental Geology under the Chinese Academy of Geological Sciences and the Natural Resources Monitoring Institute of Henan Province (Table 1).

Fig.2 Hydrogeological cross sections

2.2 Study m ethods

2.2.1 PS-InSAR

The PS-InSAR processing procedure involves several crucial steps, including image selection and registration, differential interferometric processing,of PS points selection, and inversion for surface deformation rates (Lei et al.2016; Li et al.2022).In this study, these steps were conducted as follows: (1) First, the SAR images were arranged in chronological order, and the images from December 15, 2018, were selected as the reference master images for both the northern and southern parts.These master images were used to register the remaining slave images.By setting temporal and spatial vertical baseline thresholds at 90 days and 150 days, respectively, 29 interferometric pairs were generated for each part (Fig.3); (2) Next, the topographic and flat-land errors in the interferometric images were elim inated using SRTM DEM data w ith a resolution of 90 m, form ing differential interferograms.Subsequently, stable target candidate points w ith strong scattering, known as PS points, were extracted in the study area using the amplitude dispersion index.A fter performing the phase unw rapping of PS points, the phase information representing their linear deformations was obtained and was then subtracted from the original differential interferometric phase, yielding the residual phase information; (3) Finally, the interference from nonlinear deformation and atmospheric phase was removed through filtering, and the deformation rate of a PS point along the direction of radar line-of-sight (LOS) was derived and further transformed into the vertical deformation rate using equation (1):

2.2.2 GIS spatial analysis

GIS spatial analysis utilizes geostatistics and spatial operations as its theoretical foundation, to analyze the spatial distribution and characteristics of various types of geographic data w ithin a specificarea, thereby extracting new spatial information.Commonly used analysis methods in GIS include buffer analysis, overlay analysis, equal-fan analysis, and fishnet analysis (Peng, 2008; Qu et al.2014).

Table 1 Data types, description, and sources

Fig.3 Spatial and temporal baselines of Sentinel-1A images

Based on the geological environment of the study area, as well as data from leveling measurements, groundwater regime monitoring, and the multi-source monitoring of surface deformations,this study constructed the spatial data field of surface deformations in the study area and their associated driving factors through the following steps: (1) First, GIS spatial analysis was used to conduct numerical statistics, allowing for the analysis of the spatio-temporal evolution of surface deformations; (2) Subsequently, overlay analysis and equal-fan analysis were applied to separately superimpose the driving factors of various surface deformations onto the PS-InSAR-derived calculation results.This analysis aimed to investigate the responses of surface deformations to these driving factors; (3) Finally, the indices required by the geographical detector model were gridded using fishnet tool, and the grid data on the same scale were extracted.

2.2.3 Geographical detector model

Spatially stratified heterogeneity refers to variations in a specific attribute value observed across different areas or different land use types or soil types.The geographical detector model is a statistical approach for exploring the spatial heterogeneity of geospatial feature distributions and their underlying driving mechanisms.The fundamental principle of the geographical detector model is that, when an independent variable has significant incluence on a certain dependent variable, both variables should exhibit similar spatial distributions (Wang and Xu, 2017).The geographical detector model consists of four modules, namely the factor detector, interaction detector, risk detector, and ecological detector.For this study, the factor detector and the interaction detector were primarily utilized (Wang and Xu, 2017).

(1) Factor detection

Factor detection is designed primarily to detect the spatially stratified heterogeneity of a dependent variable and quantify the extent to which various independent variables contribute to the spatial variation.It aims to determine the importance of various factors in the research context, denoted byq.A higherqvalue indicates a more substantial influence of a driving factor on the spatially stratified heterogeneity of land subsidence, with the expression as follows:

Where:his the layer of Y or X;andNdenotes the number of units in layerhand the whole area, respectively;andare the variances of value Y of layerhand the whole area,respectively, and the value range ofqis [-1,1].

(2) Interaction detection

Interaction detection can identify whether the combined effects of different driving factors enhance or weaken their individual contributions to a dependent variable, or whether the driving factors have independent influences on that variable.Specifically, theqvalues for two variables related to land subsidence rate can be calculated separately.Then, theqvalue in the case of interactions between the two variables (a new layer formed by the superimposing of layers X1and X2) can be calculated.Finally, the impact of the interactions between the two variables can be assessed by comparing theqvalues before and after considering interactions.The types of interactions and their corresponding criteria are presented in Table 2.

3 Results and discussion

3.1 InSAR monitoring results and accuracy verification

According to statistics, a total of 81,099 PS points were extracted in the study area from the images processed using the PS-InSAR technique.These PS points are distributed primarily on buildings and along the edges of road in towns and villages,suggesting relatively reliable inversion results(Fig.4).

Table 2 The interaction types and corresponding criteria

Fig.4 Distribution of leveling measurement points, PS points, and benchmark of the leveling measurement

The accuracy of the vertical deformation rates obtained through InSAR inversion was verified by comparing them with synchronous observations from 52 leveling measurement points in the study area.During data selection process, if a leveling monitoring point did not coincide with a PS point,its InSAR result was set to the average value of all PS points within a circular buffer zone with a radius of 300 m and the leveling monitoring point at its center.The errors of the two measurement methods were compared using equations (3) and(4).As shown in Fig.5, the cumulative deformation amplitude calculated using both methods had an average error of 6.61 mm, a root mean square error of 6.88 mm, and a correlation coefficient of up to 0.957, indicating a strong correlation.Based on the aforementioned error analysis results, it can be concluded that the InSAR monitoring results effectively reflect the surface deformations in the study area.

Fig.5 Verification results of InSAR monitoring

3.2 Spatio-temporal distribution of surface deformations

After applying Kriging interpolation method to process the PS-InSAR results, this study generated a map depicting the spatial distribution of surface deformation rates in the study area from January 2018 to March 2020 (Fig.6).As shown in Fig.6,the surface deformations in the study area exhibit a significantly heterogeneous spatial pattern, with deformation rates ranging from -37.1 mm/a to 8 mm/a.The land subsidence zones cover an area of 1,562.1 km2, accounting for 71.17% of the total study area.Land subsidence, as the primary form of surface deformation, is distributed primarily in the central and eastern parts of Zhongmu County and the Yellow River beach area in the northern part of this county.The surface uplift zones are concentrated in the south central part of the main urban area and the Airport Economy Zone, accounting for about 19.83% of the total study area.

Fig.6 Spatial distribution of surface deformations rate in the study area

To further analyze the spatial distribution of surface deformations in the study area, this study classified the developmental degree of surface deformations based on deformation rates and cumulative deformation amplitude according to the Specifications for Risk Assessment of Geological Hazard and the actual situation in Zhengzhou(Table 3).By setting deformation rates greater than-8 mm/a and greater than 0 mm/a as the thresholds for typical surface deformations zones, this study identified four typical land subsidence zones and two surface uplift zones, with the former consisting of A (Putian-Baisha-Liuji), B (Wantan-Dongzhang-Langchenggang), C (Dameng-Cangzhai-Guandu), and D (Huiji-Jinshui) and the latter comprising E (Gufo-Mazhai-Nancao) and F (Bagang-Zhangzhuang-Sanguanmiao).

The spatial distribution of the four typical subsidence zones is shown in Fig.6, and the details are as follows: (1) Zone A covers the main urban area and the eastern new urban area.It experiences severe land subsidence, with deformation rates ranging from -23 mm/a to -37.1 mm/a.The subsidence center lies in the groundwater source area of Baisha Town.It has the highest cumulative deformation amplitude of -88.4 mm during the monitoring period, indicating the most severe subsidence in the study area; (2) Zone B is located in the eastern part of Zhongmu County, and has its subsidence center in Cangzhai Town.Its deformation rates range from -12.8 mm/a to -26.3 mm/a, with a cumulative deformation amplitude of -57 mm during the monitoring period; (3) Zone C is locatedin the Yellow River beach area in the northern part of Zhongmu County and displays an irregular laminar pattern.It is the largest zone in the study area with moderate land subsidence ranging from-8 mm/a to -22.6 mm/a and its cumulative deformation amplitude is -33 mm during the monitoring period Its subsidence center lies in Langchenggang Town; (4) Zone D is sporadically distributed in the northern part of the main urban area, with subsidence rates ranging from -8 mm/a to -27.4 mm/a.Although this zone previously experienced the most severe land subsidence in Zhengzhou(Guan et al.2019), currently only Gangli Village in Huiji District and Shamen Village in Jinshui District face are subjected to fairly severe subsidence, with maximum cumulative deformation amplitude of -77 mm during the monitoring period.Consequently, the subsidence center of the study area has shifted from the Huiji-Jinshui area to Baisha Town in Zhongmu County.The two surface uplift zones are described as follows: (1)Zone E is located within the area where groundwater exploitation is restricted, and it falls in the water-receiving area of the South-to-North Water Diversion Project (phase II).The surface deformation rates in this zone range from +4 mm/a to+8 mm/a.Zone F is located in the northern part of Airport Economy Zone and also falls within the water-receiving area of the South-to-North Water Diversion Project (phase II).The surface deformation rates range from +1.8 nm to +5.6 mm.

Table 3 Criteria for developmental degrees of surface deformations

To analyze the temporal variations of the surface deformations in the study area, this study conducted statistical analysis on the cumulative deformation amplitude of all PS points during the monitoring period.As shown in Table 4, the maximum cumulative deformation amplitude in the land subsidence zones increased from -33.82 mm in Stage I (from January 2018 to February 2019) to-40.13 mm in Stage II (from February 2019 to March 2020), with an increase of 6.54 mm.On the other hand, the cumulative deformation amplitude in the surface uplift zones showed minimal variation between the two stages, with changes less than 1 mm, which can be considered negligible within the permissible error range of PS-InSAR.Further analysis was conducted on the zones with different developmental degrees of land subsidence (Table 5).In Stage II, the areas of zones with severe and moderate land subsidence were found to be 26.34 km2and 432.42 km2, respectively.These areas increased by 10.1 km2and 94.82 km2, respectively compared to Stage I, suggesting significant increase in amplitude.In addition, the surface uplift zones exhibited an increase in area of 33.80 km2in Stage II compared to Stage I.

As indicated by the spatial distribution and temporal variations, the surface deformations in the study area were dominated by land subsidence during the monitoring period, and the regional land subsidence exhibited a slow development over this period.

3.3 Analysis of factors driving surface deformations through GIS spatial analysis

Numerous studies have demonstrated that surface deformations are influenced by a combination of natural and anthropogenic factors.Natural factors include tectonism, lithology, and sea level rise,while anthropogenic factors comprise groundwater exploitation and urban construction (Motagh et al.2017; Wang et al.2018; Ye et al.2021).The effective stress of soils can increase or decrease depending on the groundwater level and pore water pressure, causing either compaction or rebound deformation of soils.In this case, the lithology of soil mass is a crucial factor that influences surface deformations (Motagh et al.2017).In areas with intensive urban construction and heavy traffic, the compactional deformation of stratum surfaces or underlying soil mass can be accelerated by static loads from buildings, dynamic loads induced by traffic, excavation of foundation pits, and precipitation, all of which contribute to severe land subsidence (Wang et al.2018; Zhou et al.2019).Therefore, considering both natural and anthropogenic factors, this section utilizes GIS spatial analysis to explore the relationship between surface deformations and groundwater exploitation, thickness of compressible layers, and urban construction.

Table 4 Statistics of surface deformation rates in the study area

Table 5 Statistics of changes in the area of zones with different developmental degrees of surface deformations from January 2018 to March 2020

3.3.1 Groundwater exploitation

Due to the increased proportion of diverted water supply and the implementation of restrictions on groundwater exploitation, the groundwater exploitation volume in the main urban area of Zhengzhou has shown a declining trend.Since 2017, the main urban area has achieved a positive equilibrium between groundwater extraction and recharge, leading to a gradual rise in groundwater levels in areas with restricted groundwater exploitation (Guan et al.2019; Ye, 2021).In contrast, Zhongmu County has experienced relatively severe groundwater over-exploitation and has been in a negative equilibrium between groundwater extraction and recharge since 2017, except for the years 2016 and 2018.These findings align with the spatial distribution of surface deformations.Specifically, zones with severe and moderate land subsidence are primarily located in Zhongmu County, while zones with minimal land subsidence and surface uplift are concentrated in the main urban area (Fig.7).

Fig.7 Groundwater exploitation in the study area

To further investigate the relationship between surface deformations and groundwater exploitation, this study conducted a spatial analysis using ArcGIS by overlaying the groundwater level contours of the first and second aquifer groups in 2019 onto the cumulative deformation amplitude observed during the monitoring period.The analysis revealed that the zones experiencing severe land subsidence exhibited a spatial distribution that roughly corresponded to the subsidence funnels of the moderately deep groundwater; however, their centers showed some deviation (Fig.8).This can be attributed to the following factors: The second aquifer group, with exploitation depths ranging from 80 m to 360 m, mainly consists of alluvium,diluvium, and lacustrine silty clay and clay.The permeability coefficient of these cohesive soils is three to five orders of magnitude lower than that of water-bearing sand layers.Consequently, the consolidation of cohesive soils after water release lags behind the change in the water head (Jia et al.2007).

In addition, a greater drop in groundwater levels can induce higher cumulative deformation amplitude in subsidence funnels in the second aquifer group.For example, Lihuqiao Village, experiencing the most severe land subsidence during the monitoring period, is located in the Baisha Town-Cangzhai Town-Guandu Town subsidence funnel.The center of this funnel had a groundwater level decline of 6.7 m, corresponding to a cumulative deformation amplitude of -90.96 mm over the monitoring period.Similarly, in the Guandu Town-Hansi Town area, the groundwater level of the moderately deep aquifer group decreased by 1.34 m and 1.23 m, resulting in cumulative deformation amplitudes of -44.47 mm and 22.88 mm,respectively.In contrast, monitoring wells in Mazhai Town, located in a zone with surface uplift, showed an increase in the groundwater levels of the first and second aquifer groups by 3.3 m and 2.6 m, respectively, corresponding to a surface uplift of 15 mm.These findings indicate a strong correlation between surface deformations and groundwater levels in the study area.

Fig.8 Elevation contours of groundwater level and cumulative deformation amplitude in the study area in 2019

3.3.2 Lithology

The unconsolidated deposits in the Quaternary and the Neogene strata play a crucial role in the surface deformations of the study area (Pan et al.2020;Fig.9), as highlighted in previous studies.The lithology of these deposits is a significant factor influencing the deformations of soil mass.It has been observed that cohesive soils exhibit primarily plastic deformations, while sandy soils tend to disply elastic deformations during the process of soil consolidation after water release.Moreover,the consolidation-induced deformation amplitude of the cohesive soils is typically two to three times greater than that of sandy soils under the same pressure (Guo et al.2017; Lei et al.2022).

Fig.9 Spatial distribution of cohesive soil thickness and cumulative deformation amplitude

To further investigate the relationships between lithology and surface deformations, this study conducted an analysis of the cohesive soil thickness across the study area using data from 92 boreholes.The obtained cohesive soil thickness map was then overlaid with the cumulative deformation amplitude data.The resulting map (Fig.9),reveals that the zones experiencing severe and moderate land subsidence are predominantly distributed in areas with cohesive soil thicknesses greater than 200 m.On the other hand, zones with weak land subsidence and surface uplift are mostly found in areas with cohesive soil thicknesses less than 180 m.This distribution can be attributed to the behavior of cohesive soils after water release.When subjected to continuous groundwater exploitation, cohesive soils undergo significant compaction and elastic deformation.Consequently, areas with thick cohesive soils experience severe and irreversible land subsidence.In contrast, areas characterized by thin cohesive soils and thick sandy soils exhibit lower deformation amplitude caused by groundwater exploitation.Furthermore, these areas primarily undergo elastic deformation, which can recover with a decrease in groundwater exploitation volume and an increase in groundwater level.

3.3.3 Urban construction

Using the reclassification tools in ArcGIS, this study extracted images showing the evolution of land use types in the study area from 2010 to 2019(Fig.10).Moreover, a GIS spatial analysis was conducted to analyze the relationship between the growth rate of construction land and the cumulative deformation amplitude (Fig.11).The analysis results reveal that urban expansion in Zhengzhou over the past decade has occurred in a radial pattern, originating from the main urban area and extending towards Zhongmu County.Specifically,there has been significant urban construction and expansion towards Zhengdong New District and Baisha Town, which is aligns with the spatial distribution of zones experiencing severe land subsidence with cumulative deformation amplitude greater than -20 mm, as well as zones with moderate land subsidence.In contrast, zones with minimal land subsidence with cumulative deformation amplitude less than -10 mm and zones with surface uplift are mostly distributed in the south central part of the main urban area, where the growth rate of construction land has been relatively low.

Fig.10 Land use types in the study area during 2010-2019 (1: Construction land; 2: Cultivated land; 3:Water body; 4: Grassland; 5: Forest land)

Fig.11 Evolutionary relationships between construction land and zones with land subsidence based on equal-fan analysis

Table 6 shows the corresponding relationships between the growth area of construction land and the zones with different levels of land subsidence accross various equal fans.The highest growth rates of construction land area were observed in the NW, NEE, E, and SEE directions, with area increases of 115.46 km2, 123.91 km2, 98.88 km2,and 45.56 km2, respectively, over the course of 10 years.The zones experiencing land subsidence in these four directions covered an area of about 1,154 km2, which accounted for 64% of the total subsidence area within the study area.Among these subsidence zones, areas with cumulative deformation amplitudes greater than -30 mm and -20 mm to -30 mm covered 59% and 64% of the total subsidence area in the study area, respectively.These results indicate that the significant increase in highdensity building loads and frequent engineering construction activities have accelerated local land subsidence (Zhao et al.2020).

3.4 Analysis of factors driving surface deformations using the geographical detector model

As demonstrated by spatial correlation analysis in Section 3.3, there is a clear correspondence between the cumulative deformation amplitude and individual factors such as groundwater exploitation, lithology, and urban construction.However,there is a lack of quantitative analysis of the contribution of each driving factor to the spatially stratified heterogeneity of surface deformations.In this section, the quantifiable driving factors were extracted first.Subsequently, the contributive degrees of these factors to the cumulative deformation amplitude were quantitatively described using the factor detector in the geographical detector model, which helps identify the dominant driving factors.Furthermore, the contributive degrees of various driving factors after pairwise interactions were analyzed using the interaction detector in the geographical detector model.Finally, the driving factors influencing the spatial heterogeneity of surface deformations were explored by combining the analytical results of the two detectors.

3.4.1 Selection of quantifiable driving factors

Based on the analytical results from Section 3.3,this study identified several driving factors that contribute to the spatial stratified heterogeneity of land subsidence.These factors include building density (X1), road densities (X2), groundwater levels in the first and second aquifer groups (X3and X4, respectively), and cohesive soil thicknesses in the first and second aquifer groups (X5and X6, respectively).These factors jointly constituted independent variable X (Table 7), and the cumulative deformation amplitude was taken as dependent variable Y in the model.To facilitate the analysis, the study area was divided into fishnet grids measuring 300 m × 300 m using the Fishnet tool in ArcGIS.Subsequently, the independent variable X (Table 7) and dependent variable Y were classified using the natural breaks method,and the classification results were exported as corresponding grid files (Fig.12).In total, data from 27,352 uniformly distributed fishnet grids were collected, and the X and Y values were input into the geographical detector model for further analysis.

3.4.2 Contributive degrees of driving factors

The results obtained in Section 3.2 reveal thatsurface deformations in the study area exhibit significant spatially stratified heterogeneity.To analyze the contributive degrees of various driving factors to the heterogeneity, the factor detector was employed.Table 8, illustrates the interpretation degrees (represented bypvalues) of the driving factor, where lowerpvalues indicate higher interpretation degrees.According to Table 8, all the selected driving factors exhibited significant contributive degrees (p＜ 0.05) to the spatially stratified heterogeneity of surface deformations, confirming the appropriateness of the chosen factors.The contributive degrees (q) of these driving followed the order: X4＞ X3＞ X2＞ X5＞ X6＞ X1.Among these factors, X4demonstrated the highest contributive degree with aqvalue of 0.5382, making it the most influential factor in shaping the surface deformations in the study area.X3ranked the second with aqvalue of 0.2674.X2had aqvalue of 0.0983, slightly higher than theqvalues of X5and X6, while X1exhibited the lowest contributive degree with aqvalue of 0.0067.Therefore,groundwater exploitation emerges as the primary cause of the significant spatially heterogeneity of surface deformations within the study area currently.

Table 6 Statistics of the increased area of construction land and the area of zones with different cumulative deformation amplitudes from January 2018 to March 2020

Table 7 Driving factors selected for the geographical detector model

3.4.3 Interactions of indicators

Fig.12 Spatial distribution of driving factors

Table 8 Detection results of driving factors

Based on the results obtained from the factor detector, this study analyzed the changes in the contributive degrees of various driving factors by considering their interactions using the interaction detector.The findings, as shown in Fig.13, indicate that the contributive degrees of these driving factors to the spatially stratified heterogeneity were significantly enhanced after pairwise interactions.The interaction types included two-factor enhancement and nonlinear enhancement.

As shown in Fig.13, the driving factors withqvalues greater than 0.5 after pairwise interactions included X4∩X6, X3∩X4, X2∩X4, X1∩X4, and X4∩X5.After interactions with X4, other driving factors showed the largest increased amplitudes of contributive degrees.Among these, the interaction between X4and X6demonstrate the strongest control effects on the spatial stratified heterogeneity of land subsidence in the study area, with X4∩X6having aqvalue of 0.5722.It is important to note that both X5and X6initially displayed low contributive degrees, withqvalues of 0.0396 and 0.0713, respectively.However, their contributive degrees significantly improved after interactions with X3and X4, respectively, resulting inqvalues of 0.3123 and 0.5722, respectively.

These results indicate that surface deformations in the study area are dominated by the soil deformations caused by changes in groundwater level,especially the compactional deformation of soils caused by the groundwater overexploitation from the second aquifer group.The reasons for this can be explained as follows: The second aquifer group has relatively poor recharge conditions, and the water exploited volume in this aquifer group is mostly sourced from the groundwater storage that cannot be recovered by the aquifer itself.Moreover, the second aquifer group has a large proportion of cohesive soils (Fig.2) and thus will generate larger-amplitude irreversible deformations during the continuous over-exploitation of groundwater (Guo et al.2017; Lei et al.2022).In contrast,the first aquifer group has favorable recharge conditions and can be rapidly recharged by meteoric precipitation and surface water; with the recovery of groundwater level, the elastic sandy soils in this aquifer group will rebound partially (Jia et al.2007; Lei et al.2022).Consequently, severe land subsidence occurs in areas characterized by a substantial decline in groundwater levels and relatively thick cohesive soils, such as Baisha, Dameng, and Guandu towns.In contrast, areas with significant rise in groundwater levels, low content of cohesive soils, and high content of sandy soils,such as Mazhai Town, Shifo Town, and Nancao Subdistrict, exhibit minimal land subsidence, and local surface uplift may even occur.

Fig.13 Detection results of driving factor interaction

4 Conclusion

Presently, the plain area in eastern Zhengzhou is predominantly experiencing gradual land subsidence.Moreover, the center of this regional land subsidence has shifted from the main urban area to Zhongmu County due to an imbalanced water supply structure and the rapid expansion of the construction land in the subcentral of Zhengzhou City.

The GIS spatial analysis reveals a strong spatial correlation between the cumulative deformation amplitude in the study area and factors such as groundwater exploitation, lithology, and urban construction during the monitoring period.

The spatially stratified heterogeneity of surface deformations in the study area is influenced by multiple driving factors, including building density(X1), road density (X2), the groundwater levels of the first and second aquifer groups (X3and X4,respectively), and the cohesive soil thicknesses in the first and second aquifer groups (X5and X6,respectively).Among these factors, X4is the dominant driving factor.Furthermore, the interactions between X4and X6contribute the most to the spatially stratified heterogeneity, indicating that changes in groundwater levels plays a crucial role in the soil deformation leading to the surface deformations in the study area.

Acknowledgements

This paper was supported by the China Geological Survey Project (Grant No.DD20189262; Grant No.DD20211309) and Basic Research Operations Project of the Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences (SK202206).