Detection of Land Subsidence Associated with Land Creation and Rapid Urbanization in the Chinese Loess Plateau Using Time Series InSAR : A Case Study of Lanzhou New District

Lanzhou New District is the first and largest national-level new district in the Loess Plateau region of China. Large-scale land creation and rapid utilization of the land surface for construction has induced various magnitudes of land subsidence in the region, which is posing an increasing threat to the built environment and quality of life. In this study, the spatial and temporal evolution of surface subsidence in Lanzhou New District was assessed using Persistent Scatterer Interferometric Synthetic Aperture radar (PSInSAR) to process the ENVISAT SAR images from 2003–2010, and the Small Baseline Subset (SBAS) InSAR to process the Sentinel-1A images from 2015–2016. We found that the land subsidence exhibits distinct spatiotemporal patterns in the study region. The spatial pattern of land subsidence has evidently extended from the major urban zone to the land creation region. Significant subsidence of 0–55 mm/year was detected between 2015 and 2016 in the land creation and urbanization area where either zero or minor subsidence of 0–17.2 mm/year was recorded between 2003 and 2010. The change in the spatiotemporal pattern appears to be dominated mainly by the spatial heterogeneity of land creation and urban expansion. The spatial associations of subsidence suggest a clear geological control, in terms of the presence of compressible sedimentary deposits; however, subsidence and groundwater fluctuations are weakly correlated. We infer that the processes of land creation and rapid urban construction are responsible for determining subsidence over the region, and the local geological conditions, including lithology and the thickness of the compressible layer, control the magnitude of the subsidence process. However, anthropogenic activities, especially related to land creation, have more significant impacts on the detected subsidence than other factors. In addition, the higher collapsibility and compressibility of the loess deposits in the land creation region may be the underlying mechanism of macro-subsidence in Lanzhou New District. Our results provide a useful reference for land creation, urban planning and subsidence mitigation in the Loess Plateau region, where the large-scale process of bulldozing mountains and valley infilling to create level areas for city construction is either underway or forthcoming.


Introduction
Urbanization has become an irreversible global trend [1,2], and the limited supply of available land is an important impediment for rapid urbanization, especially in the case of cities in upland areas such as the Loess Plateau in China; examples are the cities of Lanzhou and Yan'an.In August 2012, the Chinese Government State Council approved the state-level new district in Lanzhou (Lanzhou New District, LZND), in Gansu Province, where the thickest loess deposits in the world occur [3]; it is the first and largest national-level new area in the Loess Plateau region.There is a severe shortage of available land for urban development in the region, and consequently local governments have resorted to removing the tops of many high loess mountains to fill the adjacent valleys to create flat land for construction-a process of land creation.The land creation project has contributed greatly to urban construction and the local economy, and the areas have experienced rapid and significant growth of large-scale engineering construction since 2013.However, there is a strenuous debate underway regarding the environmental impacts of land creation [4,5], with the major concern that large-scale land creation may cause severe land subsidence.This subsidence poses a major threat to urban construction with the possibility of substantial damage to the urban infrastructure, including buildings, railroads and highways, and increasing the risk of flooding and accelerated erosion along land fissures [6].
The cities in the Loess Plateau are expanding rapidly, accompanied by a significant increase in large-scale land creation and urbanization.However, until now, very little deformation monitoring activity, such as point measurements (GNSS: Global Navigation Satellite System) or large-scale high-precision monitoring work, has been conducted in this region, where wholesale bulldozing of mountains is underway to build cities such as LZND.Consequently, the rates and patterns of region-wide subsidence in LZND are unknown.Synthetic Aperture Radar Interferometry (InSAR) is an ideal methodology for exploring the spatial extension, magnitude, and temporal evolution of land subsidence and its possible causes in a typical loess-underlain city associated with large-scale land creation and rapid urbanization.The technique can measure ground surface displacements over wide areas (at both regional-and local-scale) at a millimeter-scale of accuracy for long periods of time [7,8].InSAR has been widely applied in many cities for monitoring land subsidence related to human activities such as groundwater exploitation and urbanization [9][10][11][12][13][14][15][16][17].The creation of time-series of subsidence monitoring with high temporal resolution relies mainly on advanced multi-temporal interferometry methods, such as persistent scatter interferometry (PSI), SBAS (Small Baseline Subsets) and the combined analysis of the two [18][19][20].It has been determined that PSI enables the precise characterization of linear deformation (slow or very slow moving) where the "persistent scatters" (e.g., buildings) are abundant (such as urban areas), while SBAS may be more effective in rural or mountainous areas where few persistent scatters can be detected, such as variations in landcover and vegetation coverage [10,[21][22][23].To characterize the complex deformation process across areas of rapid urbanization and large-scale land creation where the radar targets are not stable or limited, the combined use of these two techniques is preferable.
The present study uses the combination of PSInSAR and SBAS techniques to quantify the subsidence phenomenon and its spatial distribution and temporal evolution in a typical loess-covered region (LZND) where wholesale bulldozing of mountains to build cities is underway.Subsequently, InSAR-derived subsidence is analyzed in conjunction with different anthropogenic and hydrogeological factors to determine the possible causes.We first present an overview of the study area and the methodology used.Then, we present the spatiotemporal evolution of land subsidence in LZND from 2003 to 2010 and from 2015 to 2016, processed by the PSInSAR technique using ENVISAT SAR images and the SBAS methodology using Sentinel-1A images.Finally, we discuss the possible causes of the observed subsidence to better understand its relationship with environmental and anthropogenic variables.In addition, the possible underlying mechanics of land subsidence are discussed.

Study Area
LZND is in the Qingwangchuan Basin in northern Lanzhou City, 38 km from the downtown area (Figure 1).It is also centrally located in relation to the three provincial capital cities of Lanzhou (Gansu Province), Xining (Qinghai Province) and Yinchuan (Ningxia Province).The region is an important portal for the opening and development of Gansu Province, and is also a significant connection between the Silk Road economic zone and the Eurasian "Continental Bridge".

Study Area
LZND is in the Qingwangchuan Basin in northern Lanzhou City, 38 km from the downtown area (Figure 1).It is also centrally located in relation to the three provincial capital cities of Lanzhou (Gansu Province), Xining (Qinghai Province) and Yinchuan (Ningxia Province).The region is an important portal for the opening and development of Gansu Province, and is also a significant connection between the Silk Road economic zone and the Eurasian "Continental Bridge".The study area is a mountain basin, the northern part of which consists of low-elevation mountains, and the remainder of low loess hills with a relative elevation of 40-60 m.The strata in the basin are mainly Quaternary loose sediments and Neogene mudstone.The composition is Quaternary alluvial silt, silty clay, sand, and gravel, underlain by Neogene mudstone and sandstone.The mountains around the basin are mainly covered by Late Pleistocene loose aeolian loess; these mountains with loose loess are the major source of material to fill the valleys to create flat land for construction.
The study area has a typical semi-arid continental monsoon climate.The average annual temperature and precipitation are 4.1 °C and 300-350 mm, respectively.The spatial and temporal distribution of precipitation in the LZND is uneven; with increasing elevation, the rainfall increases from southeast to the northwest, and is concentrated during July-September.The region has high evaporation with an average annual amount of 1880 mm.The study area is a mountain basin, the northern part of which consists of low-elevation mountains, and the remainder of low loess hills with a relative elevation of 40-60 m.The strata in the basin are mainly Quaternary loose sediments and Neogene mudstone.The composition is Quaternary alluvial silt, silty clay, sand, and gravel, underlain by Neogene mudstone and sandstone.The mountains around the basin are mainly covered by Late Pleistocene loose aeolian loess; these mountains with loose loess are the major source of material to fill the valleys to create flat land for construction.
The study area has a typical semi-arid continental monsoon climate.The average annual temperature and precipitation are 4.1 • C and 300-350 mm, respectively.The spatial and temporal distribution of precipitation in the LZND is uneven; with increasing elevation, the rainfall increases from southeast to the northwest, and is concentrated during July-September.The region has high evaporation with an average annual amount of 1880 mm.
Groundwater in Qinwangchuan is mainly distributed in the ancient trench on both sides of the basin, while in other parts of the basin the groundwater is unevenly distributed.The groundwater in the study region can be divided into two categories: The first aquifer is phreatic water in the Quaternary loose strata, and the second is confined water within the Neogene interbedded clastic rocks.The groundwater depth decreases from south to north, from more than 30 m in the ancient trench in northern Daoshuitang and Zhoujialiang, to 10-20 m in Fangjiapo-Xicao in the southern basin (Figure 1).A large-scale hydraulic project has been implemented-the "Lead into the Qin Engineering"-which consists of the large-scale transfer of water from Datong River through canals in Qinghai Province into the Qinwangchuan Basin, with an annual water diversion volume of 443 million m 3 .The irrigation area has increased significantly and the groundwater level has gradually increased [24].The project also provides the major water source for urban consumption in the LZND.Large-scale land creation projects, involving cutting of the loess mountains to infill the valleys to create flat land for construction, commenced in 2013 (Figure 2) and are concentrated mainly in the southeastern and northeastern parts of the basin (Figure 1).Due to the bareness of the aeolian loess surface in the land creation region, the potential stable scatterers were limited, which may lead to a low density of PS in the Envisat results.Groundwater in Qinwangchuan is mainly distributed in the ancient trench on both sides of the basin, while in other parts of the basin the groundwater is unevenly distributed.The groundwater in the study region can be divided into two categories: The first aquifer is phreatic water in the Quaternary loose strata, and the second is confined water within the Neogene interbedded clastic rocks.The groundwater depth decreases from south to north, from more than 30 m in the ancient trench in northern Daoshuitang and Zhoujialiang, to 10-20 m in Fangjiapo-Xicao in the southern basin (Figure 1).A large-scale hydraulic project has been implemented-the "Lead into the Qin Engineering"-which consists of the large-scale transfer of water from Datong River through canals in Qinghai Province into the Qinwangchuan Basin, with an annual water diversion volume of 443 million m 3 .The irrigation area has increased significantly and the groundwater level has gradually increased [24].The project also provides the major water source for urban consumption in the LZND.Large-scale land creation projects, involving cutting of the loess mountains to infill the valleys to create flat land for construction, commenced in 2013 (Figure 2) and are concentrated mainly in the southeastern and northeastern parts of the basin (Figure 1).Due to the bareness of the aeolian loess surface in the land creation region, the potential stable scatterers were limited, which may lead to a low density of PS in the Envisat results.

Data and Methods
To determine the spatiotemporal changes in land subsidence in the LZND, we utilized two sets of C-band SAR images, including 40 descending SLC scenes acquired by ENVISAT ASAR from June 2003 to September 2010, and 15 ascending SLC scenes from Interferometric Wide Swath (IW) mode of Sentinel-1A acquired from May 2015 to July 2016, which cover the LZND and the surrounding area (Figure 1).SRTM (Shuttle Radar Topography Mission) with 90-m resolution was used to correct the topography-related phase and geocode interferograms.To explore the spatial association between subsidence and anthropogenic factors, multi-temporal high-resolution images (16 June 2003, 11 June 2006, and 3 March 2013 Google Earth images and 17 October 2016 GF-1 satellite image) were used to interpret the land use types.In addition, data from 43 boreholes were used to interpolate the soil depth and groundwater level in the basin.In addition, a consolidation experiment was conducted on 48 and 16 soil samples from the infilled area and original hilly area, respectively, to determine soil geotechnical parameters.
In this research, the PSI and SBAS were applied to process the 40 ENVISAT ASAR images and 15 Sentinel-1A images (Table 1), respectively, both incorporated into the ENVI platform through the SARSCAPE 5.2 module (http://www.sarmap.ch/).The reasons we used different InSAR time-series methods to process the images for different time periods are as follows: (1) from 2003 to 2010, the PSInSAR was selected because of the large number of images, no or slow deformation, and a stable radar signal (e.g., from buildings) which did not vary substantially across the region because the urbanization process was relatively slow during this period; and (2) from 2015 to 2016, potential permanent scatters such as buildings changed progressively due to rapid urbanization, and few

Data and Methods
To determine the spatiotemporal changes in land subsidence in the LZND, we utilized two sets of C-band SAR images, including 40 descending SLC scenes acquired by ENVISAT ASAR from June 2003 to September 2010, and 15 ascending SLC scenes from Interferometric Wide Swath (IW) mode of Sentinel-1A acquired from May 2015 to July 2016, which cover the LZND and the surrounding area (Figure 1).SRTM (Shuttle Radar Topography Mission) with 90-m resolution was used to correct the topography-related phase and geocode interferograms.To explore the spatial association between subsidence and anthropogenic factors, multi-temporal high-resolution images (16 June 2003, 11 June 2006, and 3 March 2013 Google Earth images and 17 October 2016 GF-1 satellite image) were used to interpret the land use types.In addition, data from 43 boreholes were used to interpolate the soil depth and groundwater level in the basin.In addition, a consolidation experiment was conducted on 48 and 16 soil samples from the infilled area and original hilly area, respectively, to determine soil geotechnical parameters.
In this research, the PSI and SBAS were applied to process the 40 ENVISAT ASAR images and 15 Sentinel-1A images (Table 1), respectively, both incorporated into the ENVI platform through the SARSCAPE 5.2 module (http://www.sarmap.ch/).The reasons we used different InSAR time-series methods to process the images for different time periods are as follows: (1) from 2003 to 2010, the PSInSAR was selected because of the large number of images, no or slow deformation, and a stable radar signal (e.g., from buildings) which did not vary substantially across the region because the urbanization process was relatively slow during this period; and (2) from 2015 to 2016, potential permanent scatters such as buildings changed progressively due to rapid urbanization, and few persistent scatters could be detected in the land creation region located in rural or mountainous areas.In addition, the SBAS reduces decorrelation in areas with high deformation rates.Many studies have been conducted to validate the consistency between PSI and SBAS results, and they demonstrate the utility of combining both methods to detect subsidence and landslide movements [25][26][27].Therefore, it is appropriate to apply both methods to detect and compare ground deformation in LZND.

PSInSAR
PSInSAR identifies the permanent scatterers by analyzing the amplitude and phase variation in a series of interferograms.These scatterers dominate the returning phase in their pixels which are brighter than those of background scatterers [28].These stable scatterers can be anthropogenic (e.g., buildings and metallic structures) or natural objects (e.g., rock outcrops) [29].A detailed description of the PSInSAR technique is given in Ferretti et al. [7].The applied PSI processing steps are summarized as follows: master image selection, co-registration, differential interferogram computation, PS candidates (PSC) selection, estimation of displacement and topographic correction for PSC, atmospheric phase screen (APS) estimation and removal, PS selection, mean displacement rate estimation and topographic errors removal, PS points displacement time series analysis and average deformation estimation [30,31].

SBAS
Although the minimum number of interferograms to properly estimate the deformation signal for PSInSAR is usually 15 [32], the accuracy could increase with the increasing number of images [23].To obtain more reliable and denser results in the land creation region where the subsidence rates were relatively high, the Small Baseline Subset method was used to process the 15 scene Sentinel-1A images.The SBAS method relies on processing the interferogram with a shorter baseline and fewer decorrelations to determine the surface displacement through time [10].For SBAS processing, four main steps were involved.
The first key step of the SBAS algorithm is the generation of a set of multilook differential interferograms, which begins with a selection of SAR image pairs with spatial and temporal constraints (smaller than 210 m and less than 180 days) that are represented in a "connection graph" to form interferograms (Figure 3).Sixty-two differential interferograms were generated after the removal of the topographic phase.Multi-looking with 5 × 1 in range and azimuth direction and the Goldstein filtering method [33] was applied to improve the signal/noise ratio of the interferograms.The next step is to conduct phase unwrapping of the wrapped phase signal calculated from the selected interferograms.The phase of the coherent pixels was unwrapped by performing the Minimum Cost Flow approach (MCF) [34,35].Several interferogram pairs with low coherence and unwrapped phase error were discarded [36].The residual phase content and phase ramps were calculated to correct the unwrapped phase by selecting and refining stable Ground Control Points (GCPs).In this research, we selected GCPs based on the following criteria: (1) temporal coherence value over 0.75; (2) the location was far away from areas with large displacement based on interferogram interpretation and field survey; and (3) deformation close to zero confirmed by interpretation of unwrapped phase and field survey [37][38][39].Consequently, we used 51 GCPs in both orbital refinement and displacement conversion steps.Then, the preliminary displacements were estimated using a linear model and residual topography was also removed.The Singular Value Decomposition (SVD) method was applied to search the least squares solution for each coherent pixel as well as to estimate the nonlinear deformation.Finally, we derived the time-series deformation after subtracting the estimated atmospheric artifacts and orbital ramps.

SBAS
Although the minimum number of interferograms to properly estimate the deformation signal for PSInSAR is usually 15 [32], the accuracy could increase with the increasing number of images [23].To obtain more reliable and denser results in the land creation region where the subsidence rates were relatively high, the Small Baseline Subset method was used to process the 15 scene Sentinel-1A images.The SBAS method relies on processing the interferogram with a shorter baseline and fewer decorrelations to determine the surface displacement through time [10].For SBAS processing, four main steps were involved.
The first key step of the SBAS algorithm is the generation of a set of multilook differential interferograms, which begins with a selection of SAR image pairs with spatial and temporal constraints (smaller than 210 m and less than 180 days) that are represented in a "connection graph" to form interferograms (Figure 3).Sixty-two differential interferograms were generated after the removal of the topographic phase.Multi-looking with 5 × 1 in range and azimuth direction and the Goldstein filtering method [33] was applied to improve the signal/noise ratio of the interferograms.The next step is to conduct phase unwrapping of the wrapped phase signal calculated from the selected interferograms.The phase of the coherent pixels was unwrapped by performing the Minimum Cost Flow approach (MCF) [34,35].Several interferogram pairs with low coherence and unwrapped phase error were discarded [36].The residual phase content and phase ramps were calculated to correct the unwrapped phase by selecting and refining stable Ground Control Points (GCPs).In this research, we selected GCPs based on the following criteria: (1) temporal coherence value over 0.75; (2) the location was far away from areas with large displacement based on interferogram interpretation and field survey; and (3) deformation close to zero confirmed by interpretation of unwrapped phase and field survey [37][38][39].Consequently, we used 51 GCPs in both orbital refinement and displacement conversion steps.Then, the preliminary displacements were estimated using a linear model and residual topography was also removed.The Singular Value Decomposition (SVD) method was applied to search the least squares solution for each coherent pixel as well as to estimate the nonlinear deformation.Finally, we derived the time-series deformation after subtracting the estimated atmospheric artifacts and orbital ramps.

Results
In the present study, 48 pairs of interferograms images were generated for the PSI analysis (Figure 3a), and a coherence coefficient of 0.75 was used to maintain confidence in the generation of high quality displacements.For the SBAS calculation, two thresholds of 210 m and 180 days were used for the spatial and temporal baseline, respectively (Figure 3b), to derive surface displacement rates in line-of-sight (LOS) across the LZND.The PSI-and SBAS-derived deformation along the LOS (dLOS) was transformed to the vertical plane (dV) using the following relationship: 40] based on the ENVISAT and Sentinel-1A incidence angle, respectively.Thus, a positive value

Results
In the present study, 48 pairs of interferograms images were generated for the PSI analysis (Figure 3a), and a coherence coefficient of 0.75 was used to maintain confidence in the generation of high quality displacements.For the SBAS calculation, two thresholds of 210 m and 180 days were used for the spatial and temporal baseline, respectively (Figure 3b), to derive surface displacement rates in line-of-sight (LOS) across the LZND.The PSI-and SBAS-derived deformation along the LOS (d LOS ) was transformed to the vertical plane (d V ) using the following relationship: [9,40] based on the ENVISAT and Sentinel-1A incidence angle, respectively.Thus, a positive value represents an uplift motion and the negative motion rates' sign indicates movement away from the satellite (i.e., land subsidence).The results indicate that surface deformation in the LZND is evident, primarily in the range of −17.2 mm/year to 13.3 mm/year for 2003-2010, and in the range of −55 mm/year to 21 mm/year for 2015-2016.The precision of the estimated displacement rates was determined based on the standard deviation of the deformation velocity [23]; it ranged 0-2 mm/year and 0-8 mm/year for 2003-2010 and 2015-2016, respectively, and confirms the reliability of the estimated displacements.There are 180,782 PS points and 1,028,684 Coherent Targets (CTs) over the study region, with an area of approximately 36 × 25 km 2 , consisting of 201 PS points km −2 and 1143 CTs km −2 .The much higher spatial density of CTs may be attributed to the fact that high quality interferograms were guaranteed due to the small baseline strategy that was introduced in interferometric generation.The reliability and accuracy of the deformation results were further validated through field investigation since there was no ground-based leveling data for this region.

Spatial Distribution of Land Subsidence
Deformation maps from 2003 to 2010 and from 2015 to 2016 in the LZND were generated to illustrate the spatial distribution of land subsidence (Figure 4).It can be seen in Figure 4 that the spatial distribution of the subsidence in the LZND changed substantially.From 2003 to 2010, there are three small areas with zero or very low subsidence rates, ranging 0-17.2 mm/year (Figure 4a).The major subsidence was concentrated in the central-eastern part of Zhongchuan Town, northwestern Xicha Town (Figure 4a), and the central urban area of Qingchuan Town (Figure 4a).From 2015 to 2016, the extent of land subsidence expanded substantially, and the three areas with the highest localized subsidence rates were Sidun-Fangjiapo district (Figure 4b), Shengli-Liudun district (Figure 4b) and Baojiayao-Yuchuan district (Figure 4b).It can be seen in Figure 4 that during 2003-2010 the two subsidence areas (Figure 4a) in the northwestern part of Xicha Town and the eastern part of Zhongchuan Town apparently expanded to form a large area of subsidence (Figure 4b).Within the area, this resulted in uneven settlement, and the subsidence rate also substantially increased, with maximum subsidence rates for these three major subsidence centers of 55.0 mm/year, 32.4 mm/year and 22.6 mm/year, respectively.The subsidence was mainly related to earth movement and valley infilling work, as well as the rapid increase in construction.It is noticeable that the subsidence is also spreading to the southeastern part of Qingchuan Town due to land creation.In addition, substantial land subsidence has commenced in the Baojiayao-Yuchuan district, which was generally stable during 2003-2010, whereas a maximum subsidence rate of 22.6 mm/year was observed for 2015-2016.According to the spatial distribution of land subsidence, the three major land subsidence areas are Sidun-Fangjiapo district, Shengli-Liudun district and Baojiayao-Yuchuan district, which are analyzed in detail.

Land Subsidence in Sidun-Fangjiapo District
This district is in the main urban area of the LZND (Figure 4b).From 2003 to 2010, the main settlement area was in Sidun and Fangjiapo, and other locations; subsidence in the main urban area was not obvious, with zero or very low subsidence rates.The average and maximum subsidence rates were 2.40 mm/year and 4.06 mm/year, respectively.Between 2015 and 2016, the two subsidence areas of Sidun and Fangjiapo were combined into a large area of subsidence.The average subsidence rate increased substantially to 10.0 mm/year, and the maximum subsidence rate was 55.0 mm/year.The land subsidence in this area was mainly related to land creation, and the external load resulted from the construction of high-rise buildings and road construction in the main urban area, which has formed a series of unevenly distributed subsidence features, such as ground fissures, building cracks and foundation subsidence.

Land Subsidence in Shengli-Liudun District
This subsidence area is mainly a recently-developed industrial region (Figure 4b) and the spatial expansion of land subsidence is essentially consistent with engineering activity.From 2003 to 2010, the subsidence rate was relatively low, with an average rate of 2.8 mm/year and a maximum rate of 17.2 mm/year.Subsidence was mainly concentrated in the center of Qinchuan Town, and thus the spatial extent is small.However, the average subsidence rate for 2015-2016 was 10.4 mm/year.There was an apparent subsidence center in the land creation region of Wudun village, with a maximum subsidence rate of 32.4 mm/year.The spatial distribution of subsidence was extended substantially compared with 2003-2010, resulting from large-scale engineering activity such as land creation, economic and technological development and road construction.Subsidence clearly expanded to the eastern loess hills (e.g., Liudun village and Wudun village) and the western region (Shengli Village).

Land Subsidence in Baojiayao-Yuchuan District
This region consists mainly of suburban areas (Figure 4b).There was minimal land subsidence in the area from 2003 to 2010, with an average subsidence rate of only 0.78 mm/year, and a maximum rate of 1.98 mm/year.Between 2015 and 2016, there was evident land subsidence in the area, with an average subsidence rate of 9.9 mm/year and a maximum rate of 22.6 mm/year.The aeolian Malan loess, with abundant macropores and strong collapsibility, is widely distributed in the area and is the major source for valley infilling.The loose structure of the loess may be the main cause of land subsidence, and engineering activities, including earth moving and infilling, and road construction, were probably the major anthropogenic factors responsible for subsidence.

Land Subsidence in Shengli-Liudun District
This subsidence area is mainly a recently-developed industrial region (Figure 4b) and the spatial expansion of land subsidence is essentially consistent with engineering activity.From 2003 to 2010, the subsidence rate was relatively low, with an average rate of 2.8 mm/year and a maximum rate of 17.2 mm/year.Subsidence was mainly concentrated in the center of Qinchuan Town, and thus the spatial extent is small.However, the average subsidence rate for 2015-2016 was 10.4 mm/year.There was an apparent subsidence center in the land creation region of Wudun village, with a maximum subsidence rate of 32.4 mm/year.The spatial distribution of subsidence was extended substantially compared with 2003-2010, resulting from large-scale engineering activity such as land creation, economic and technological development and road construction.Subsidence clearly expanded to the eastern loess hills (e.g., Liudun village and Wudun village) and the western region (Shengli Village).

Land Subsidence in Baojiayao-Yuchuan District
This region consists mainly of suburban areas (Figure 4b).There was minimal land subsidence in the area from 2003 to 2010, with an average subsidence rate of only 0.78 mm/year, and a maximum rate of 1.98 mm/year.Between 2015 and 2016, there was evident land subsidence in the area, with an average subsidence rate of 9.9 mm/year and a maximum rate of 22.6 mm/year.The aeolian Malan loess, with abundant macropores and strong collapsibility, is widely distributed in the area and is the major source for valley infilling.The loose structure of the loess may be the main cause of land subsidence, and engineering activities, including earth moving and infilling, and road construction, were probably the major anthropogenic factors responsible for subsidence.

Evolution of Land Subsidence from 2003 to 2016
As can be seen in Figure 4, there was a substantial change in maximum subsidence rates during the study period, and the evolution of spatial patterns is obvious.In general, the central and south-eastern parts of the urban area are experiencing substantially increasing land subsidence rates.In addition, new land subsidence commenced in the western and eastern land creation region during 2015-2016.The calculated deformation shows that substantial land subsidence is spreading to the west and east where a large-scale land creation project has been conducted, with subsidence rates up to 55.0 mm/year.
We focused on the three sectors of land creation, building construction and irrigation, which are likely implicated in land subsidence.The deformational behavior of these specific sectors is especially important for understanding the process of evolution of subsidence in LZND.These sectors are designated the "Southern sector" (southeast of Zhongchuan Town), "Central sector" (central urban zone in Zhongchuan Town) and "Western sector" (west of Qingchuan Town).
The Southern sector (Figure 5) is in the center of Fangjiapo land creation district which was the major earthmoving and infilling region.It experienced an average deformation rate of −1.18 mm/year during 2003-2010.From 2015 to 2016, the subsidence rate increased substantially with an average rate of 35.17 mm/year.The thickness of the compressible deposits in this area exceeds 25 m.The displacement time series exhibits continuous deformation activity, and it is obvious that the deformation rate increased substantially.The results are consistent with a field investigation which revealed numerous ground fissures, sinkholes and wall cracks in the infilling area (Figure 5d-f).

Evolution of Land Subsidence from 2003 to 2016
As can be seen in Figure 4, there was a substantial change in maximum subsidence rates during the study period, and the evolution of spatial patterns is obvious.In general, the central and south-eastern parts of the urban area are experiencing substantially increasing land subsidence rates.In addition, new land subsidence commenced in the western and eastern land creation region during 2015-2016.The calculated deformation shows that substantial land subsidence is spreading to the west and east where a large-scale land creation project has been conducted, with subsidence rates up to 55.0 mm/year.
We focused on the three sectors of land creation, building construction and irrigation, which are likely implicated in land subsidence.The deformational behavior of these specific sectors is especially important for understanding the process of evolution of subsidence in LZND.These sectors are designated the "Southern sector" (southeast of Zhongchuan Town), "Central sector" (central urban zone in Zhongchuan Town) and "Western sector" (west of Qingchuan Town).
The Southern sector (Figure 5) is in the center of Fangjiapo land creation district which was the major earthmoving and infilling region.It experienced an average deformation rate of −1.18 mm/year during 2003-2010.From 2015 to 2016, the subsidence rate increased substantially with an average rate of 35.17 mm/year.The thickness of the compressible deposits in this area exceeds 25 m.The displacement time series exhibits continuous deformation activity, and it is obvious that the deformation rate increased substantially.The results are consistent with a field investigation which revealed numerous ground fissures, sinkholes and wall cracks in the infilling area (Figure 5d-f).The Central sector (Figure 6), located in the main urban region of Zhongchuan Town, which is experiencing major urbanization, was generally stable during 2003-2010, but subsidence apparently occurred in 2015-2016, with an average rate of −10.37 mm/year.From 2003 to 2010, this area was residential with low-rise buildings.Numerous roads began to be constructed across the town, and many high-rise buildings have appeared since 2013.The traffic density and associated vibration, as well as loading from construction, could induce substantial unevenly-distributed subsidence (Figure 6e-g).
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 23 The Central sector (Figure 6), located in the main urban region of Zhongchuan Town, which is experiencing major urbanization, was generally stable during 2003-2010, but subsidence apparently occurred in 2015-2016, with an average rate of −10.37 mm/year.From 2003 to 2010, this area was residential with low-rise buildings.Numerous roads began to be constructed across the town, and many high-rise buildings have appeared since 2013.The traffic density and associated vibration, as well as loading from construction, could induce substantial unevenly-distributed subsidence (Figure 6e-g).The Western sector, in Qinchuan Town, was generally stable and there was minimal subsidence during 2003-2010, with an average deformation rate of 0.31 mm/year.During 2015-2016, the subsidence rate increased substantially to 8.55 mm/year (Figure 7a-c).During the field survey, we noted numerous ground fissures in areas of cultivated land (Figure 7d), subsidence of the adjacent housing foundations (Figure 7e) and wall cracking (Figure 7f).This has formed a hazard chain: irrigation leads to cultivated land settlement and subsequently housing foundation subsidence and eventually the development of cracks in housing.Although residential properties are concentrated in this area, it is relatively far from the main urban region and the degree of urbanization is relatively limited.Comparison of remote sensing images spanning several years indicated that this region comprises mainly low-rise buildings and there was no obvious change in the urban structure pattern between 2013 and 2016, and that the built area is adjacent to agricultural land.We infer that the increase in irrigation and local wastewater drainage may be the main causes of the uneven subsidence pattern in this area.The Western sector, in Qinchuan Town, was generally stable and there was minimal subsidence during 2003-2010, with an average deformation rate of 0.31 mm/year.During 2015-2016, the subsidence rate increased substantially to 8.55 mm/year (Figure 7a-c).During the field survey, we noted numerous ground fissures in areas of cultivated land (Figure 7d), subsidence of the adjacent housing foundations (Figure 7e) and wall cracking (Figure 7f).This has formed a hazard chain: irrigation leads to cultivated land settlement and subsequently housing foundation subsidence and eventually the development of cracks in housing.Although residential properties are concentrated in this area, it is relatively far from the main urban region and the degree of urbanization is relatively limited.Comparison of remote sensing images spanning several years indicated that this region comprises mainly low-rise buildings and there was no obvious change in the urban structure pattern between 2013 and 2016, and that the built area is adjacent to agricultural land.We infer that the increase in irrigation and local wastewater drainage may be the main causes of the uneven subsidence pattern in this area.

Effect of Large-Scale Land Creation and Construction Infilling
Large-scale land creation, consisting of the removal of the tops of mountains to infill valleys, commenced in 2013 due to the construction of LZND.During 2013-2016, the land creation region expanded rapidly and covered an area about 70 km 2 .We analyzed the relationship between the land creation area, construction infilling (mainly in the urban zone) and subsidence.It can be seen in Figure 8 that the spatial distribution of subsidence is consistent with the land creation and construction infilling areas: 61.5% of the total CTs with displacement rates below −8 mm/year, indicating subsidence, are in these areas.However, a greater number of CTs, with displacement rates less than −8 mm/year, indicating substantial subsidence, occurred in the land creation region (36.7%)compared with the construction infilling area (24.8%).Moreover, the subsidence rate in the land creation region was substantially higher than in the construction infilling or other areas; in addition, a localized subsidence rate over 21 mm/year was only detected in the former (Figure 8a,b).This is especially apparent in the major land creation areas, such as Yuchuan, Liudun and Fangjiapo (Figure 8).

Effect of Large-Scale Land Creation and Construction Infilling
Large-scale land creation, consisting of the removal of the tops of mountains to infill valleys, commenced in 2013 due to the construction of LZND.During 2013-2016, the land creation region expanded rapidly and covered an area about 70 km 2 .We analyzed the relationship between the land creation area, construction infilling (mainly in the urban zone) and subsidence.It can be seen in Figure 8 that the spatial distribution of subsidence is consistent with the land creation and construction infilling areas: 61.5% of the total CTs with displacement rates below −8 mm/year, indicating subsidence, are in these areas.However, a greater number of CTs, with displacement rates less than −8 mm/year, indicating substantial subsidence, occurred in the land creation region (36.7%)compared with the construction infilling area (24.8%).Moreover, the subsidence rate in the land creation region was substantially higher than in the construction infilling or other areas; in addition, a localized subsidence rate over 21 mm/year was only detected in the former (Figure 8a,b).This is especially apparent in the major land creation areas, such as Yuchuan, Liudun and Fangjiapo (Figure 8).The maximum subsidence rates during 2015-2016 in Liudun village and Fangjiapo were 32.4 mm/year and 55.0 mm/year, respectively.Land creation in the hilly loess region, especially associated with valley infilling work, may potentially modify the physical properties (void ratio, water content and structure) and other characteristics (compressibility, collapsibility) of the naturally-compressible loess (Table 2), due to remolding during the filling process (See detailed discussion at Section 5.3).Therefore, in the infilling area, the soil, with a loose structure and in a primary consolidated state, may experience substantial self-induced compression.In addition, wetting (rainfall or irrigation) and external loading (engineering vehicles and high-rise buildings) could also accelerate the consolidation and collapse of the infilled soil, which could eventually lead to substantial land subsidence.We infer that the land creation is the predominant triggering factor determining the spatiotemporal pattern of land subsidence in LZND, and thus we suggest that the construction work should not be extended to the infilled areas until sufficient time has lapsed for ground stabilization to occur through natural settling.
The maximum subsidence rates during 2015-2016 in Liudun village and Fangjiapo were 32.4 mm/year and 55.0 mm/year, respectively.Land creation in the hilly loess region, especially associated with valley infilling work, may potentially modify the physical properties (void ratio, water content and structure) and other characteristics (compressibility, collapsibility) of the naturally-compressible loess (Table 2), due to remolding during the filling process (See detailed discussion at Section 5.3).Therefore, in the infilling area, the soil, with a loose structure and in a primary consolidated state, may experience substantial self-induced compression.In addition, wetting (rainfall or irrigation) and external loading (engineering vehicles and high-rise buildings) could also accelerate the consolidation and collapse of the infilled soil, which could eventually lead to substantial land subsidence.We infer that the land creation is the predominant triggering factor determining the spatiotemporal pattern of land subsidence in LZND, and thus we suggest that the construction work should not be extended to the infilled areas until sufficient time has lapsed for ground stabilization to occur through natural settling.Urbanization, including the magnitude of the imposed building load, has been demonstrated to be another factor involved in promoting consolidation and high displacement rates [8].In recent years, with the rapid increase in construction in LZND, especially with the movement of a zone of economic and technological development into the region, large-scale and intensive high-rise building construction has increased dramatically (Figure 9a-c).There were 290 building groups, covering an area of ~42.77 km 2 before 2003, while the number of buildings increased rapidly between 2003 and 2016 due to urbanization.In particular, during 2013-2016, the construction area expanded rapidly, increasing to 419 building groups and covering 24.65 km 2 in area, which was twice the rate of increase during 2003-2013.The rapid increase in the number of buildings has resulted in a significant load being imposed on the ground within a short interval of time.Since the soils in LZND are typically soft sediments with high compressibility, the rapid increase in ground loading may impose increased stresses on this weak stratum, which could induce compression of macropores, eventually resulting in obvious land subsidence.The distribution of low-rise and high-rise buildings in LZND from 2003 to 2016 (Figure 9a) shows that the extent of the built area was obviously expanding.It expanded rapidly from the central urban area to the surrounding area, and the land subsidence exhibits a high correlation with the distribution of building land, and the area with a high building density (the center of Zhongchuan Town) corresponds to a relatively high subsidence rate (Figure 9c).From 2003 to 2010, the areas of localized subsidence (even with small subsidence rates) were concentrated in the central urban zone where the building density was high.In addition, from 2015 to 2016, 13.5% of the total CTs with displacement rates below −8 mm/year, indicating subsidence, were in the building construction area, whereas about one-half of the CTs indicating subsidence associated with building construction correspond with the buildings constructed during 2013-2016.The consolidation process caused by the construction load can be divided into two phases [41]: a primary consolidation stage, in which the greatest subsidence is recorded; and a secondary consolidation stage, where only cohesive and organic soil can produce a much smaller amount of compression with slow creep [8].These two consolidation stages are highly time dependent and are largely determined by the thickness, hydraulic conductivity, drainage conditions and compressibility of the soft soil layer, as well as the duration and magnitude of the building load applied.A typical example of this process is building groups "B" and "C" in the center of Zhongchuan Town and in the industrial district, respectively, where the rate of subsidence was relatively high during 10 March 2015-10 March 2016.This is typical of rapid deformation in the primary consolidation stage.The mean subsidence rate time series shows a decrease in the rate during 10 March 2016-7 March 2016, which may reflect a reduction in the rate of consolidation involved in the transition from the primary to the secondary phase.However, the times series deformation of building groups "A" and "D" during 10 May 2016-10 July 2016 shows an increased subsidence rate, which may be related to a continuous primary consolidation process.The subsidence rate of building groups "A" and "D" exceeds 9 mm/year and can be ascribed to the geotechnical characteristics of the loose loess with high compressibility and low hydraulic conductivity.In addition, the dimensions of the structure, which comprise a large area of high-rise buildings, could transfer a high loading to the ground, inducing rapid subsidence.It can be concluded that in LZND the subsidence rate increases with the increase in the speed and scale of construction, as well as the building density.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 23 that in LZND the subsidence rate increases with the increase in the speed and scale of construction, as well as the building density.

Effects of the Rapid Increase in Road Construction
During the construction of the transport network in LZND, with the occurrence of abrupt changes in ground elevation along the road routes, many subgrades have been filled with various thicknesses of loose Quaternary loess.Under the action of dynamic loads such as cyclic traffic loading and the consequent vibrations, as well as pile driving, consolidation of the soft soil has occurred, which could significantly exacerbate the progression of land subsidence.It can be seen in Figure 10 that the extent and density of the road network increased rapidly from 2013 to 2016.In particular, in

Effects of the Rapid Increase in Road Construction
During the construction of the transport network in LZND, with the occurrence of abrupt changes in ground elevation along the road routes, many subgrades have been filled with various thicknesses of loose Quaternary loess.Under the action of dynamic loads such as cyclic traffic loading and the consequent vibrations, as well as pile driving, consolidation of the soft soil has occurred, which could significantly exacerbate the progression of land subsidence.It can be seen in Figure 10 that the extent and density of the road network increased rapidly from 2013 to 2016.In particular, in the southern and eastern part of LZND, where the main urbanization zone (Zhongchuan Town) and the land creation region are located, a dense road network was constructed between 2013 and 2016.The density of road networks within LZND was 0.92 km/km 2 and 1.90 km/km 2 before 2013 and during 2013-2016, respectively.It can be seen in Figure 10 that the distribution of the dense road network is well correlated with the distribution of subsidence.Many examples of unevenly distributed subsidence were detected by time-series InSAR along the road networks; in addition, the number of CTs with displacement rates exceeding the threshold (−8 mm/year) indicate that subsidence decreased with increasing distance to the road network, and higher rates of subsidence were registered at locations closer to the road network.All of the subsidence was within 800 m of the road networks, and 62.4% of the total CTs with displacement rates below −8 mm/year (indicating subsidence) are concentrated within 200 m of the road network.In addition, subsidence at rates over 21 mm/year was only observed within 400 m of the road (see Figure 10).A series of ground subsidence occurrences is distributed on road pavements or embankments, which was consistent with ground fissures, road pavement cracks and embankment settlement noted during the field investigation.subsidence decreased with increasing distance to the road network, and higher rates of subsidence were registered at locations closer to the road network.All of the subsidence was within 800 m of the road networks, and 62.4% of the total CTs with displacement rates below −8 mm/year (indicating subsidence) are concentrated within 200 m of the road network.In addition, subsidence at rates over 21 mm/year was only observed within 400 m of the road (see Figure 10).A series of ground subsidence occurrences is distributed on road pavements or embankments, which was consistent with ground fissures, road pavement cracks and embankment settlement noted during the field investigation.

Relationship between Land Subsidence and Geology
Surface geological characteristics, such as lithology, are an important factor influencing subsidence [10].As can be seen in Figure 11, there is a very clear lithological control on the spatial distribution of subsidence in LZND.The Holocene alluvial gravel soil, loess-like silt and late Pleistocene aeolian loess are the most widely distributed deposits in LZND, which covered 47.9% and 42.5% of the area in this region, respectively.Almost all of the subsidence (66.4% and 33.4% of the total CTs with displacement rates below −8 mm/year, indicating subsidence for the Holocene and late Pleistocene deposits, respectively) is distributed in the areas of these two types of soft deposit with a loose structure.However, no subsidence is observed in areas with lithologies of sand and mudstone, conglomerate, phyllite and slate.Furthermore, the higher rates of subsidence are mainly concentrated in the aeolian loess-covered region where extensive efforts have been made to create a flat ground surface by levelling the hills and infilling the valleys; in particular, more CTs (474 vs. 381) with subsidence rates over 15 mm/year were located in the late Pleistocene aeolian loess areas than in areas with Holocene deposits.It should be noted there was either zero subsidence or minor subsidence during 2003-2010, even though the same geological conditions occurred; and that it was the large-scale land creation and rapid urbanization that were responsible for the evident land subsidence in areas underlain by these soft deposits during 2013-2016.

Relationship between Land Subsidence and Compressible Deposits
The distribution of land subsidence is also affected by the presence and thickness of compressible deposits [14].Quaternary loose sediments are the main compressible deposits in the Qingwangchuan Basin.Differences in the thickness of the layer of Quaternary sediments is an important geological factor determining the occurrence and development of uneven land subsidence in the LZND.To analyze the spatial relationship between the thickness of compressible deposits and subsidence, we overlaid the map of land subsidence for 2015-2016 on that of the thickness of compressible deposits (Figure 12a).It is evident that the thickness of compressible deposits controls the spatial distribution of land subsidence.Generally, the greater the thickness of the soft soil layer the higher subsidence rate (Figure 12b-d), and this indicates that subsidence is relatively large under the action of external loading and self-compression, since the Quaternary loose sediments are highly compressible.As can be seen in Figure 12a, land subsidence has occurred mainly (96.3% of the total CTs with displacement rates below −8 mm/year) in areas where the thickness of compressible deposits ranges from around 15-45 m; while the maximum subsidence rate is observed in the southeastern basin, including Fangjiapo and Zhaojiapu, with thicker compressible deposits.In contrast, in the center and northern parts of the basin, with shallow compressible deposits (from 0-15 m), the subsidence rate was either much lower or non-existent, which was 3.7% of the total CTs with displacement rates below −8 mm/year occurring in these areas.It should be noted that obvious subsidence was also observed in Shengli-Liudun where a relatively shallow soft soil is distributed; this may be ascribed to the continuous and high-intensity land creation work being conducted in the area.

Relationship between Land Subsidence and Groundwater
It is generally recognized that groundwater withdrawal was the major driver of subsidence in several major cities, such as Mexico City, Beijing and Shanghai [9,11,[42][43][44].There has been no large-scale groundwater exploitation in LZND in recent years, and the groundwater level has increased due to the "Lead into the Qin Engineering" hydraulic diversion project [24].As can be seen in Figure 13, there is a weak correlation between subsidence and the spatial distribution of groundwater: a portion of the settlement area corresponds with a shallow groundwater depth, such as Fangjiapo and Shengli-Liudun; however, high rates of subsidence are also observed in western Qingchuan and Xicha where the groundwater depth was relatively deep.It can be concluded that groundwater fluctuations did not play a significant role in land subsidence in LZND.However, the subsequent increase in the phreatic water level in the Quaternary loose sediments, due to the "Lead into the Qin Engineering" hydraulic diversion project, may have caused structural collapse of the shallow loess, and this could have resulted in subsidence [42].Conversely, in the near-future, continuous rapid urbanization could significantly increase the demand for water, and large-scale groundwater exploitation, which lowers the groundwater level and may also induce subsidence.Therefore, InSAR combined with groundwater monitoring should be conducted to support sustainable urbanization in the LZND.

Relationship between Land Subsidence and Groundwater
It is generally recognized that groundwater withdrawal was the major driver of subsidence in several major cities, such as Mexico City, Beijing and Shanghai [9,11,[42][43][44].There has been no large-scale groundwater exploitation in LZND in recent years, and the groundwater level has increased due to the "Lead into the Qin Engineering" hydraulic diversion project [24].As can be seen in Figure 13, there is a weak correlation between subsidence and the spatial distribution of groundwater: a portion of the settlement area corresponds with a shallow groundwater depth, such as Fangjiapo and Shengli-Liudun; however, high rates of subsidence are also observed in western Qingchuan and Xicha where the groundwater depth was relatively deep.It can be concluded that groundwater fluctuations did not play a significant role in land subsidence in LZND.However, the subsequent increase in the phreatic water level in the Quaternary loose sediments, due to the "Lead into the Qin Engineering" hydraulic diversion project, may have caused structural collapse of the shallow loess, and this could have resulted in subsidence [42].Conversely, in the near-future, continuous rapid urbanization could significantly increase the demand for water, and large-scale groundwater exploitation, which lowers the groundwater level and may also induce subsidence.Therefore, InSAR combined with groundwater monitoring should be conducted to support sustainable urbanization in the LZND.

Effects of Soil Collapsibility and Compressibility
To better investigate the effect of land creation on subsidence in the LZND, the geotechnical characteristics of samples from the area of land infilling and the original loess hilly area were examined (Table 2).The results show that the loess fill has a very similar grain-size distribution to the original loess, which is typical silty clay.In addition, the dry density of the infilled loess is very close to that of the original loess, which indicates that the translocated soil was substantially compacted during the process of land creation.However, the average void ratio of the infilled loess is much larger than that of the original loess, indicating that variations in the microstructure probably occurred during the infilling process.Therefore, the infilling loess is more likely to be self-consolidating under the influence of gravity or the external load.The consolidation test shows that the infilling loess and the original loess are both strongly collapsible soils and exhibit self-induced collapsibility.However, the infilling loess has a higher compressibility than the original loess, with coefficients of compressibility of 0.45 and 0.22, respectively.Therefore, under the external dynamic and static load, the infilling loess is more prone to settling.It can be inferred that the repeated compaction of the soil layer during the infilling process was insufficient to eliminate completely the danger of the collapse of the soft loess deposits.Conversely, the compressibility and collapsibility of the infilling loess increased compared with the original loess; therefore, additional material, such as sodium silicate, lime and fly ash piles, which could potentially reduce the collapse potential of the loess, may be an optional stabilization method during the process of large scale infilling [45,46].

Conclusions
We have investigated and analyzed land subsidence in the LZND, a region of rapid land creation and urban expansion.Subsidence monitoring was conducted using two different time series InSAR techniques for processing the radar data: PSInSAR analysis was applied to ENVISAT datasets for 2003-2010 (prior to large-scale construction), and SBAS was applied to the Sentinel 1A datasets for 2015-2016 (during large-scale construction).The spatial distribution and temporal evolution of land subsidence are identified and mapped and we have analyzed the relationships between multiple anthropogenic and hydrogeological factors and land subsidence.In addition, the possible effects of collapsibility and compressibility of the soil on subsidence are discussed.The main findings are as follows: (1) The spatiotemporal pattern of land subsidence has changed substantially in LZND.There was no obvious subsidence between 2003 and 2010, but, from 2015 to 2016, the extent of land subsidence in the urban area evidently expanded to the region of earth movement and infilling.The subsidence rate increased substantially and the maximum subsidence was 55 mm/year from 2015-2016.
(2) The spatial expansion and rate of increase of subsidence support the conclusion that large-scale earth movement and infilling and associated urban construction and rapid urbanization are the major anthropogenic drivers controlling the spatiotemporal evolution of subsidence in LZND.The highest subsidence rate was recorded in the infilling area.In addition, the associated increase in the density of the transport and irrigation network also exacerbated the subsidence process.Large-scale land creation was the predominant factor which determined the spatiotemporal pattern of land subsidence in this region, and thus construction work should not be extended to infilled regions without the elapse of a sufficient interval for the ground to stabilize.
There is a clear geological control (i.e., related to lithology and the compressibility of the substrate) on the spatial distribution of subsidence.All the subsiding areas detected are in areas of Holocene alluvial gravel soil, loess-like silt and late Pleistocene aeolian loess.We observed a correspondence between the subsidence rate and the thickness of the soft soil layer, with a higher subsidence rate evident in the thicker, more compressible deposits.It should be emphasized that it is major anthropogenic activity, such as large-scale land creation and rapid urbanization, that was responsible for causing obvious land subsidence in areas of soft deposits between 2013 and 2016.Moreover, groundwater fluctuations did not play a significant role in land subsidence in LZND, since the groundwater level has increased in recent years due to a large-scale hydraulic project.
(3) The infilling loess has a much greater degree of collapsibility and compressibility than the original loess, which is the possible mechanical cause of the greater subsidence in the areas of land creation.
Our observations have significant implications for effective land creation, urbanization planning and land subsidence management decisions in loess regions where mountains are bulldozed to increase the area available for urban development.In the LZND, following urban expansion, a series of large-scale construction projects has been initiated which are resulting in serious subsidence.Thus, the integration of InSAR with traditional monitoring techniques and consolidation modeling is necessary to implement effective land creation schemes and facilitate sustainable urban construction at a regional scale.Furthermore, stabilization agents should be added to and mixed with the infilling loess during the process of large-scale valley infilling work to improve its performance.

Figure 1 .
Figure 1.Location of Lanzhou New District (LZND) as covered by the Satellite image obtained on 14 October 2015.Blue and red rectangles indicate the areas covered by the Envisat and Sentinel-1A images, respectively.

Figure 1 .
Figure 1.Location of Lanzhou New District (LZND) as covered by the Satellite image obtained on 14 October 2015.Blue and red rectangles indicate the areas covered by the Envisat and Sentinel-1A images, respectively.

Figure 2 .
Figure 2. Land creation in LZND: (a) infilling work in the loess hilly region; and (b) repeated dynamic compaction during the infilling process.

Figure 2 .
Figure 2. Land creation in LZND: (a) infilling work in the loess hilly region; and (b) repeated dynamic compaction during the infilling process.

Figure 3 .
Figure 3. Spatial-temporal baselines of the generated interferograms: (a) Envisat PSI analysis; and (b) Sentinel-1A SBAS analysis.Each dot represents an image, and each line represents an interferogram.

Figure 3 .
Figure 3. Spatial-temporal baselines of the generated interferograms: (a) Envisat PSI analysis; and (b) Sentinel-1A SBAS analysis.Each dot represents an image, and each line represents an interferogram.

Figure 4 .
Figure 4. InSAR-derived average velocity map in the vertical direction over the LZND: (a) linear deformation trend of Envisat for 2003-2010 using PSI analysis, where the PS points were superimposed on the Google Earth image obtained on 3 March 2013; and (b) vertical velocity map of Sentinel-1A for 2015-2016 using the SBAS technique, where the CTs were superimposed on the GF-1 satellite image acquired on 17 October 2016.

Figure 4 .
Figure 4. InSAR-derived average velocity map in the vertical direction over the LZND: (a) linear deformation trend of Envisat for 2003-2010 using PSI analysis, where the PS points were superimposed on the Google Earth image obtained on 3 March 2013; and (b) vertical velocity map of Sentinel-1A for 2015-2016 using the SBAS technique, where the CTs were superimposed on the GF-1 satellite image acquired on 17 October 2016.

Figure 5 .
Figure 5.Typical subsidence observed in the land creation area located in the Southern sector: (a) Average 2003-2010 velocity map; the PSs were superimposed on the Google Earth image obtained on 3 March 2013.(b) Average 2015-2016 velocity map; the CTs were superimposed on the GF-1 satellite image acquired on 17 October 2016.Red rectangles in (a,b) indicate the location of the PSs and Coherent Targets selected, respectively.(c) Vertical displacement time-series and corresponding linear rates of subsidence.Envisat (blue) and Sentinel-1A (red) time series are shown.The gap between 2011 and 2014 was filled by a prediction trendline (grey dashed line in (c)).Images of: (d) ground fissures and sinkhole; (e) wall cracks; and (f) road settlement.

Figure 5 .
Figure 5.Typical subsidence observed in the land creation area located in the Southern sector: (a) Average 2003-2010 velocity map; the PSs were superimposed on the Google Earth image obtained on 3 March 2013.(b) Average 2015-2016 velocity map; the CTs were superimposed on the GF-1 satellite image acquired on 17 October 2016.Red rectangles in (a,b) indicate the location of the PSs and Coherent Targets selected, respectively.(c) Vertical displacement time-series and corresponding linear rates of subsidence.Envisat (blue) and Sentinel-1A (red) time series are shown.The gap between 2011 and 2014 was filled by a prediction trendline (grey dashed line in (c)).Images of: (d) ground fissures and sinkhole; (e) wall cracks; and (f) road settlement.

Figure 6 .
Figure 6.Typical subsidence observed in the main urban area located in the Central sector: (a) Average 2003-2010 velocity map; the PS points were superimposed on the Google Earth image obtained on 3 March 2013.(b) Average 2015-2016 velocity map; the CTs were superimposed on the GF-1 satellite image acquired on 17 October 2016.Red stars in (a,b) indicate the location of the PSs and CTs selected, respectively.(c) Vertical displacement time-series and corresponding linear rates of subsidence.Envisat (blue) and Sentinel-1A (red) time series are shown.The gap between 2011 and 2014 was filled by a prediction trendline (grey dashed line in (c,d)).Images of: (d) fence tilt; (e) sidewalk subsidence; and (f) foundation cracks.

Figure 6 .
Figure 6.Typical subsidence observed in the main urban area located in the Central sector: (a) Average 2003-2010 velocity map; the PS points were superimposed on the Google Earth image obtained on 3 March 2013.(b) Average 2015-2016 velocity map; the CTs were superimposed on the GF-1 satellite image acquired on 17 October 2016.Red stars in (a,b) indicate the location of the PSs and CTs selected, respectively.(c,d) Vertical displacement time-series and corresponding linear rates of subsidence.Envisat (blue) and Sentinel-1A (red) time series are shown.The gap between 2011 and 2014 was filled by a prediction trendline (grey dashed line in (c,d)).Images of: (e) fence tilt; (f) sidewalk subsidence; and (g) foundation cracks.

23 Figure 7 .
Figure 7.Typical subsidence observed in the main irrigation area located in the Western sector which has formed a geohazard chain: (a) Average 2003-2010 velocity map; the PSs were superimposed on the Google Earth image obtained on 3 March 2013.(b) Average 2015-2016 velocity map; the Coherent Targets were superimposed on the GF-1 satellite image acquired on 17 October 2016.Red stars in (a,b) indicate the location of the PSs and CTs selected, respectively.(c) Vertical displacement time-series and corresponding linear rates of subsidence.Envisat (blue) and Sentinel-1A (red) time series are shown.The gap between 2011 and 2014 was filled by a prediction trendline (grey dashed line in (c)).Images of: (d) irrigation-induced ground fissures; (e) foundation settlement due to cultivated land subsidence; and (f) wall cracks resulted from foundation subsidence.

Figure 7 .
Figure 7.Typical subsidence observed in the main irrigation area located in the Western sector which has formed a geohazard chain: (a) Average 2003-2010 velocity map; the PSs were superimposed on the Google Earth image obtained on 3 March 2013.(b) Average 2015-2016 velocity map; the Coherent Targets were superimposed on the GF-1 satellite image acquired on 17 October 2016.Red stars in (a,b) indicate the location of the PSs and CTs selected, respectively.(c) Vertical displacement time-series and corresponding linear rates of subsidence.Envisat (blue) and Sentinel-1A (red) time series are shown.The gap between 2011 and 2014 was filled by a prediction trendline (grey dashed line in (c)).Images of: (d) irrigation-induced ground fissures; (e) foundation settlement due to cultivated land subsidence; and (f) wall cracks resulted from foundation subsidence.

Figure 8 .
Figure 8. Land creation and construction infilling in relation to subsidence in LZND: (a) subsidence rate during 2015-2016 and different land reclamation types; (b) subsidence rate and cumulative number of CTs for different land reclamation types; (c-f) enlarged views for the three major land creation zones; (c,e,g) images extracted from the Google Earth satellite obtained on 3 March 2013; and (b,f,h) images extracted from the GF-1 satellite acquired on 17 October 2016.

Figure 8 .
Figure 8. Land creation and construction infilling in relation to subsidence in LZND: (a) subsidence rate during 2015-2016 and different land reclamation types; (b) subsidence rate and cumulative number of CTs for different land reclamation types; (c-f) enlarged views for the three major land creation zones; (c,e,g) images extracted from the Google Earth satellite obtained on 3 March 2013; and (b,f,h) images extracted from the GF-1 satellite acquired on 17 October 2016.

Figure 9 .
Figure 9. Localized subsidence in the main urban zone in relation to building construction: (a) Average subsidence rate from 2015-2016 and buildings constructed at different periods.The CTs are overlain on a 2016 GF-1 satellite image.(b,c) Enlarged views for the central Zhongchuan Town and newly-constructed industrial region of Xicha Town, respectively.(d) Envisat and Sentinel-1A time series for buildings A, B, C and D.

Figure 9 .
Figure 9. Localized subsidence in the main urban zone in relation to building construction: (a) Average subsidence rate from 2015-2016 and buildings constructed at different periods.The CTs are overlain on a 2016 GF-1 satellite image.(b,c) Enlarged views for the central Zhongchuan Town and newly-constructed industrial region of Xicha Town, respectively.(d) Envisat and Sentinel-1A time series for buildings A, B, C and D.
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 23 the southern and eastern part of LZND, where the main urbanization zone (Zhongchuan Town) and the land creation region are located, a dense road network was constructed between 2013 and 2016.The density of road networks within LZND was 0.92 km/km 2 and 1.90 km/km 2 before 2013 and during 2013-2016, respectively.It can be seen in Figure 10 that the distribution of the dense road network is well correlated with the distribution of subsidence.Many examples of unevenly distributed subsidence were detected by time-series InSAR along the road networks; in addition, the number of CTs with displacement rates exceeding the threshold (−8 mm/year) indicate that

Figure 10 .
Figure 10.Localized subsidence in relation to road networks in LZND.Figure 10.Localized subsidence in relation to road networks in LZND.

Figure 10 .
Figure 10.Localized subsidence in relation to road networks in LZND.Figure 10.Localized subsidence in relation to road networks in LZND.

Figure 11 .
Figure 11.InSAR-derived average annual subsidence rates in LZND, with surface lithology units from the 1:100,000 geological map superimposed.

Figure 11 .
Figure 11.InSAR-derived average annual subsidence rates in LZND, with surface lithology units from the 1:100,000 geological map superimposed.

Figure 12 .
Figure 12.Relationship between the thickness of compressible deposits with subsidence in LZND: (a) isoclines of compressible deposits thickness and average subsidence rate between 2015 and 2016; (b) geological cross section and subsidence rate profile of A-A1; (c) geological cross section and subsidence rate profile of B-B1; and (d) geological cross section and subsidence rate profile of C-C1.Light grey rectangles in (b-d) indicate that the subsidence rate increased with the increasing thickness of compressible deposits.

Figure 12 .
Figure 12.Relationship between the thickness of compressible deposits with subsidence in LZND: (a) isoclines of compressible deposits thickness and average subsidence rate between 2015 and 2016; (b) geological cross section and subsidence rate profile of A-A1; (c) geological cross section and subsidence rate profile of B-B1; and (d) geological cross section and subsidence rate profile of C-C1.Light grey rectangles in (b-d) indicate that the subsidence rate increased with the increasing thickness of compressible deposits.

Figure 13 .
Figure 13.Distribution of land subsidence rate and groundwater contour map in LZND.Figure 13.Distribution of land subsidence rate and groundwater contour map in LZND.

Figure 13 .
Figure 13.Distribution of land subsidence rate and groundwater contour map in LZND.Figure 13.Distribution of land subsidence rate and groundwater contour map in LZND.

Table 1 .
Envisat and Sentinel data for Lanzhou New District.

Table 2 .
Main physical-mechanical properties of soil test samples.N denotes the number of samples on which measurements were made in each experiment or soil test, and ± values indicate the SD from the mean.