Spatial and temporal subsidence characteristics in Wuhan city (China) during 2015-2019 inferred from Sentinel-1 SAR Interferometry

. Ground subsidence is regarded as one of the most common geohazards accompanied with the rapid urban expansion in recent years. In the last two decades, Wuhan located in the alluvial Jianghan Plain has experienced great urban expansion with increased subsidence issues, i.e. soft foundation subsidence and karst collapse. Here we investigated subsidence rates in 15 Wuhan city with 2015-2019 Sentinel-1 SAR images. We found that the overall subsidence over Wuhan region is significantly correlated with the distribution of engineering geological subregions (EGSs). We further validated the InSAR measurements with better than 5 mm accuracy by comparing with levelling measurements. Subsidence centres in Qingling-jiangdi, Houhu, Qingshan and Dongxihu area were identified with displacement rates of approximately 30 mm/yr. Our results demonstrated that the dominant driven factor is ongoing constructions and the subsidence centres shifted with construction intensities. 20 Qingling-Jiangdi area in our study is a well-known site of karst collapse. We find the nonlinear subsidence


Introduction
Nowadays, land subsidence has become a serious problem along with the rapid urban expansion (Xue et al. 2005).Land subsidence events have been reported in major cities all over the world (e.g.Shanghai (Perissin et al. 2012), Beijing (Hu et al. 2019, Zhou et al. 2019), Seville (Ruiz-Constá n et al. 2017), Texas (Kim et al. 2019), Hanoi (Dang et al., 2014) and Jakarta (Ng et al. 2012, Chaussard et al. 2013)).Over 50 cities in China have been suffering from land subsidence due to various factors (Yin et al. 2005).The major causing factors of urban land subsidence are extensive pumping of groundwater (Xue et al. 2005, Yin et al. 2005), ground fissures (Zhao et al. 2018) and tectonic faults (Xue et al. 2005, Hu et al. 2019) which threat the normal operation of urban systems and people's daily lives.The approximately accumulated economic losses caused by land subsidence reached up to 450-500 billion RMB during 1949~ 2005in China (Yin et al. 2005).Therefore, great efforts need to be made to monitor and reduce the risks induced by land subsidence and related issues.
As the largest city in central China, Wuhan has experienced significant urban expansion in the last two decades (Tan et al. 2014).The ongoing constructions of high rise buildings and metro-lines all over the metropolitan area have induced serious subsidence (Zhou et al. 2017).As a result, the subsidence areas in Wuhan have remarkably extended.Recent years, more than 40 communities in Jiang'an District, Jianghan District, Qiaokou District and Wuchang District in Wuhan City have experienced ground subsidence, resulting in cracks in buildings and municipal roads in the varying degrees (Guan et al. 2016).Although over 300 benchmarks have been set to understand the subsidence trends (Zhou et al. 2017), the coverage of the benchmarks are too sparse to capture the global picture of deformation patterns.Thus, the time-series synthetic aperture radar interferometry (InSAR) techniques (e.g.Bai et al. 2016, Costantini et al. 2016, Jiang et al., 2021) making use of the stable pixels in synthetic aperture radar (SAR) images can make up for this limitation.High resolution TerraSAR-X and COSMO-Skymed SAR images are used to investigate the spatial and temporal subsidence of Wuhan urban area during 2009-2010 (Bai et al. 2016), 2013-2014(Costantini et al. 2016), 2013-2015(Bai et al. 2019), 2012-2019(Jiang et al., 2021).They found the subsidence velocity in Houhu area reached over -70mm/yr and are mainly correlated with construction activities on quaternary soft clay and carbonate rocks areas (Costantini et al. 2016, Bai et al. 2019).Similar results were also obtained by Envisat ASAR (2008)(2009)(2010), ALOS PALSAR (2007-2010), and Sentinel-1 (2015-2019) Small baseline subset (SBAS) InSAR analysis (Han et al., 2020), Radarsat-2 SBAS InSAR analysis from 2015-2018 (Zhang et al. 2019) and medium resolution Sentinel-1 images persistent scatterer interferometry (PSI) (Benattou et al. 2018) or SBAS InSAR (Zhou et al. 2017) analysis from 2015~2017.
According to geological investigations (Guan et al. 2016, Li et al. 2019), the subsidence might be correlated with engineering geological zones which are seldom studied.
In this study, 113 Sentinle-1 SAR images from April 2015 to September 2019 covering the Wuhan metropolitan area are analysed with SBAS-InSAR method to investigate the distributions of potential risks.Comparison between InSAR and levelling measurements are conducted to validate the reliability of our measurements.We found that spatial subsidence patterns are correlated with distributions of engineering geological zones in the first terrace in Wuhan city.Relationships between time series subsidence and rainfall/ river level are also discussed.Our study area covers the Wuhan metropolitan area as indicated by the dashed rectangle in Fig. 1(a).The terrain is low and flat with maximum elevation of 240 m.Over 95% of our study area are covered by Quaternary layers (Deng et al. 1991, Xu 2016, Li et al. 2019) with diverse lithology, including gravel, sand (coarse sand, fine sand, silt), sub-sand, sub-clay, clay, muck etc. (Table S1).Moreover, the area of burial or covered dissolution carbonate rocks (a.k.a Karst) is 1091.51km 2 , accounting for 12.85% of the total area of Wuhan (Tu et al. 2019, Zheng et al. 2019)  anthropogenic activities and natural forces (Tu et al. 2019).Twenty-seven of the karst collapses occurred after 2005 and only two of them are caused by natural factors (Tu et al. 2019).Geological investigations are conducted to understand the engineering conditions for urban construction or subsidence mechanisms (Li et al. 2019, Zheng et al. 2019)

Datasets
The  3 Time-series Sentinel-1 SAR interferometry analysis

Interferomeric processing
Due to the significant variations of Doppler centroid frequency in the burst, high level co-registration accuracy with better than 0.001 pixels are required for the interferomeric analysis (Jiang, 2020).In our study, the image obtained at 26 th , November 2017 was selected as the primary image.A GPU accelerated interferometric SAR processor (Yu et al., 2019) was used to process the Sentinel-1 dataset.An AW3D DSM and Sentinels Precise Orbit Determination (POD) service were first used for geometric co-registration between consecutive bursts.Then, a network-based enhanced spectral diversity approach was employed to estimate time series azimuth shifts (Jiang, 2020).Burst de-ramping, re-ramping and resample were then carried out to resample bursts.Individual bursts can then be merged into seamless SLCs.Images with temporal baseline less than 60 days are combined to generate 368 differential interforgrams as shown in Fig. 2.

Time series displacement retrieval
The SBAS-InSAR makes use of point-like targets which remain high level of coherence over a long temporal period or slow decorrelated pixels which will remain coherent in a short period (Hooper 2008, Jiang andGuarnieri 2020).Amplitude dispersion value was used to initially select the candidates that can be used to extract useful signals (Ferretti et al. 2001).Phase stability analysis was firstly performed on these candidates.The final pixels used for the displacement rate estimation were determined by temporal coherence threshold with 0.3.Then, 3D phase unwrapping were further performed to retrieve continuous phase in the spatial and temporal dimension (Hooper and Zebker 2007).Generally, there are orbital phase ramps, 120 DEM residual, atmospheric phase and deformation phase signals in the unwrapped phase.In this paper, the phase ramps were estimated with a bilinear model (Shi et al. 2016).The DEM residual phase was estimated by the linear relationship between topographic error and baseline.Temporally high-pass and spatially low-pass filters were employed to remove the atmospheric phase.Once sources of errors are mitigated, the subsidence rate and time-series displacement history were retrieved by a leastsquares adjustment.We converted the LOS displacement to the vertical direction by means of dividing by the cosine of the 125 incidence angle as Zhou et al. (2017) did in previous study.The workflow of our study was given in Fig. 3.

Mean subsidence velocity map
Figure 4 shows the subsidence rate derived from the Sentinel-1 SAR data.Given the dense man-made objects with stable backscattering signals over long time periods over urban scenarios, a total of 8,628,652 coherent pixels are selected.The maximum displacement rate is about 40 mm/yr.There are many factors (e.g.city constructions, karst landforms and soft soils) that cause ground subsidence in Wuhan.We can notice that the positions of subsidence in Wuhan metropolitan area are mainly distributed within EGSs that are composed of compressible soft soils in the first and second terrace shown in Fig. 1(b).We find the new localized subsidence centres have emerged and old subsidence centres identified by previous studies have stabilized with decreased displacement rate.For example, the uplift areas in the right bank of Yangtze River and the subsidence centre at Jianghan and Qiaokou district caused by construction of metro lines 1 and 2 during 2009-2010 (Bai et al. 2016) have almost stabilized in this study.Localized subsidence centres in Jiang'an, Qingshan, Hongshan, Wuchang, Hanyang and Dongxihu are identified which mostly cause by the intense anthropogenic activities, such as construction of metro lines and new buildings (Bai et al. 2019).

Houhu area
Houhu area was originally lake or lake beach with underlying layers composed of highly compactable muck or muck soil with thickness ranging from 10~30 m (Sun et al. 2019).Due to activities such as ponds filling and construction of embankments, the surface was covered by soft artificial fill.Serious subsidence has occurred in this area due to the urban development in recent years (Wuhan Bureau of Natural Resources and Planning 2018).Long-term analysis of high resolution COSMO-Skymed SAR images indicated Houhu area experienced seriously subsidence from 2012 to 2019 with a maximum rate of 100 mm/yr (Jiang et al., 2021).Cracks shown in Fig. S1(a) are very widely identified on newly constructed buildings in Houhu area.The subsidence of Houhu area located in Jiang'an district is very subtle during 2009-2010 with subsidence rate of -5~5 mm/yr (Bai et al. 2016) while a maximum subsidence rate with 86 mm/yr centred at longitude 114.30° and latitude 30.6° was detected in 2013-2015 due to the construction of metro line 3 (Bai et al. 2019).However, the current velocity decreased to ~ 10 mm/yr during 2015-2019 indicated by HH3 and a new subsidence centre where HH1 located with displacement velocity of ~30 mm/yr was identified which coincided with previous study by Han et al. (2020) as shown in Fig. 8(a).The localized subsidence centres shift with the urbanization progress.Fig. 8(c) and (d) shows the ©Google Earth TM images acquired at April 2012 and October 2019.The white rectangles indicated the newly built-up areas with maximum subsidence velocity of 20~30 mm/yr.The accumulative subsidence of HH1~HH4 is given in Fig. 9.The subsidence trends are all nonlinear which might closely correlate with the construction activities in the surrounding areas.The most serious accumulative subsidence occurred at HH1 which exceeded 100 mm during April 2015 to September 2019.

Qingshan area
The subsidence rate and corresponding EGS map over Qingshan area are shown in Fig. 10(a) and (b).The maximum subsidence rate in this area reached over 20 mm/yr.It is clear that the subsidence is mainly distributed in the EGSs covered by muck soil or soft clay.The subsidence centre located in QS1 agrees with previous study that used RADARSAT-2 images from 2015 to 2018 (Zhang et al. 2019).The subsidence centre (QS2) detected in our study is also identified from the results with TerraSAR-X datasets during 2013 -2015 (Bai et al. 2019).Field investigations identified cracks on buildings in this area as shown in Fig. S1(b).Fig. 10(c) and (d) show the ©Google Earth TM images of QS2 acquired at July 2014 and October 2019 respectively, in which the white polygons highlight the newly constructed areas.The construction activities in the surrounding areas of QS2 might induce the subsidence which is also the same for QS1 in Fig. S2(a) and (b).The accumulative subsidence in QS1 and QS2 presents significant nonlinear subsidence trends in Fig. 11.The subsidence rate might be affected by the construction intensities.Meanwhile, we also notice that burial carbonatite EGSs can be found over Qingshan area, which are mainly composed of dolomitic limestone and argillaceous limestone with high content of magnesium carbonate and argillaceous limestone, and difficult to dissolve (Xu, 2016).The subsidence rate is around 5 mm/yr.Meanwhile, the susceptibility of karst collapse is considered to be low (Wuhan Bureau of Natural Resources and Planning 2018).

Groundwater pumping induced subsidence
Uneven settlements caused by excessive pumping of groundwater are common phenomenon in rapid developing cities (Chaussard et al. 2013, Wuhan Bureau of Natural Resources andPlanning 2018).Thus, the groundwater extractions are restricted in Wuhan since 2013 (Chen, 2016).However, dewatering process is generally carried out during construction of deep foundation pits (Li et al. 2013, Cui et al. 2018).As a result, groundwater level might gradually decline in the surrounding areas and subsidence occurs in the meantime.During 2014, there were approximately 10,473 building sites and 539 municipal infrastructure construction sites in Wuhan (Xu, 2016).As shown in Figs. 6, 8, 10 and 12, subsidence caused by construction activities is widely distributed in Wuhan.Taking QS1 as example, QS1 is about 150 m away from the Baoye Centre which was constructed since November 2015 in Fig. S2(a S1) which are soft clay or muck soil in the upper part and sandy soil or sandy gravel in the lower part.The groundwater level is highly correlated with the river levels of Yangtze River or the Han River, especially in the first terrace (Li et al. 2013, Chen 2016).Water level correlated displacement was also observed in the first terrace in previous studies using TerraSAR-X data (Bai et al. 2016).

Impact factors of karst subsidence
Natural factors (e.g.groundwater level variations and extreme rainfall) and anthropogenic activities will cause karst collapse (Wang et al. 2020, Tu et al. 2019).The groundwater level variations induced by river level or rainfall may create negative pressure in fissures near bedrock surface of the karst caves.In general, the groundwater in karst area can be divided into the perched water, pore confined water, fracture karst water and fracture water from the top to bottom layers (Wang et al. 2020).
The perched water and pore confined water are directly connected with the river level (Tu et al. 2019).In statistics, karst collapse directly induced by pumping of groundwater has not occurred since 2001 due to the strict control by the government (Xu, 2016).However, anthropogenic engineering activities, e.g.foundation engineering or drilling, cause the occurrences of karst collapse frequently (Zheng et al. 2019).Also, the karst collapses are closely correlated with the Yangtze River or the Han River water level fluctuation (Chen, 2016).

Conclusions
We obtained the subsidence rate map over Wuhan region from Sentinel-1 imageries acquired from 2015 to 2019.Our results were validated with levelling measurements with accuracy better than 5 mm, indicating that InSAR is an effective tool in monitoring ground subsidence with dense measurement points in urban areas.Our study revealed that the overall subsidence trends agree well with the distribution of EGSs covered by soft soils.The rapid urban development is the dominant impact factor of subsidence.Dewatering process of deep foundation pits and the corresponding consolidation lead to serious subsidence, e.g.Qingling-Jiangdi area in Fig. 6, Houhu area in Fig. 8 and Qingshan area in Fig. 10.The subsidence centres shift with the intensity of urban constructions.Furthermore, we found that Qingling-Jiangdi area located in the first terrace is suffering from karst subsidence with velocity from 20~30 mm/yr, which brings great threats to peoples' daily lives.
Comparison with precipitation and water level data, the karst subsidence might be triggered by extreme rainfall in July 2016 and the nonlinear characteristic subsidence is highly correlated with the rainfall with lags since then.Nowadays, the advent of the European Space Agency's Copernicus program and the upcoming NASA/ISRO SAR mission provides unprecedented opportunities for continuous radar mapping of the Earth with enhanced revisit frequency, which help understand the underlying driving factors and detect anomalies in subsidence under the background of urban sprawl.Automatic near-real time InSAR processing and anomalies detecting methods should be developed to timely identify the potential risks.

Figure 1 :
Figure 1: (a) Topography of our study area, (b) engineering geological subregions (EGSs) over Wuhan metropolitan area and distribution of historical karst collapse.Modified from Li et al., (2019) and Zheng et al., (2019) and (c) zoomed view of the yellow rectangle in Fig. 1(b) showing the locations of historical karst collapses and levelling points.
as shown in Fig. 1(b).Thirty-eight karst collapses (a.k.a sinkholes) marked with blue circles in Fig. 1(b) and (c) were recorded in Wuhan between 1931-2018 caused by consolidation time and poor water permeability etc.The thickness of soft soil in the first terrace generally ranges from 1 ~ 18.5 m and the maximum value reaches 37 m while the thicknesses of soft soils is relatively low in the other EGSs ranging from 1 ~ 15 m (Wuhan Bureau of Natural Resources and Planning 2018).
Sentinel-1 satellite constellation conducted by the European Space Agency (ESA) are composed of Sentinel-1A launched on 3 April 2014 and Sentinel-1B launched on 25 April 2016.It is the first sensor utilizing the Interferometric Wide swath (IW) Mode as the main acquisition mode characterized by 5 by 20 meters resolution with swath width of 250 km.113 ascending track Sentinel-1A IW SAR images between 11 th April 2015 and 29 th September 2019 were acquired as shown in Fig. 2. The ALOS World 3D 30 m (AW3D30) DSM released by the Japan Aerospace Exploration Agency (JAXA) (Takaku et al. 2016) are used for co-registration, differential interferogram generation and geocoding.Measurements from 38 levelling points marked with triangles in Fig. 1(b) obtained at 10 th September 2016, 10 th March 2017, 10 th October 2017 and 10 th May 2018 are used to validate our InSAR measurements.

Figure 2 :
Figure 2: Graph of the temporal network used for InSAR time-series analysis.

Figure 3 :
Figure 3: Workflow of our study.

Figure 4 :
Figure 4: Subsidence rate over Wuhan metropolitan area between 2015 and 2019.The boundaries of the first terrace EGZ and second terrace EGZ are indicated by grey lines.The background is hillshaded AW3D DSM.

Figure 5 :
Figure 5: InSAR measurements versus levelling measurements.(a) InSAR (20160827~20170309) vs levelling (20160910~20170310) (b) InSAR (20170309~ 20171009) vs levelling (20170309 ~ 20171010) (c) InSAR (20171009 ~ 20180513) vs levelling (20171010~ 20180510) To quantify the results, we further compared InSAR with levelling measurements.We divided the levelling data into 3 time periods and compared them with InSAR data measured at closest dates as shown in Fig.5.The levelling points are mainly distributed in the Wavy hillocky EGZ as shown in Fig. 1(c) with low displacement rates which agree with InSAR measurements.The mean, root mean square error (RMSE) and median of the difference between InSAR and levelling measurements indicated our InSAR results reached millimetre accuracy.The statistical metrics in Fig. 5(a) which are slightly larger than these in Figs.5(b) and (c) might be caused by longer time coverage of InSAR than levelling.

Figure 6 :
Figure 6: (a) Subsidence velocity of Jiangdi-Qingling area during 2015 -2019.The rectangle represent the location of (c) and (d).(b) is the corresponding EGS map.The legend is the same as Fig. 1(b) while the grey circles are the historical karst collapses.The distribution of metro lines are freely available from ©Open Street Map contributors 2020.Distributed under a Creative Commons BY-SA License.(c) and (d) are ©Google Earth TM images covering point QL3 acquired from August 2013 and December 2019. 175

Figure 8 :
Figure 8: (a) Subsidence velocity of Houhu area during 2015 -2019.The rectangle represent the location of (c) and (d).(b) is the corresponding EGS map.The legend is the same as Fig. 1(b).(c) and (d) are ©Google Earth TM images covering point HH1 acquired 180

Figure 10 :
Figure 10: (a) Subsidence velocity of Qingshan area.The rectangle represent the location of (c) and (d).(b) is the corresponding EGS map.The legend is the same as Fig. 1(b).(c) and (d) are ©Google Earth TM images covering point QS2 acquired at July 2014 and 220

Figure 11 :
Figure 11: Time series subsidence of QS1 and QS2 marked in Fig. 10(a).4.6 Dongxihu areaThe subsidence rate over Dongxihu area is shown in Fig.12, in which displacements is mainly distributed over the first terrace 225 lacustrine EGS and first terrace alluvial EGS with flat terrain.The maximum displacement is up to ~30 mm/yr.The subsidence of the carbonatite EGSs in this area was less than 10 mm/yr.The accumulative subsidence of D1 and D2 shows nonlinear displacement behaviours and reaches approximately 80 mm between 2015 and 2019 in Fig.13.The ongoing constructions, together with the artificial loading are responsible for subsidence over this region.Fig.12(c) and (d) show the ©Google Earth TM images covering D1 and D2 acquired at August 2017 and December 2017.We roughly labelled the ground changes over this 230 period with white polygons.We can see that three land cover conversions occur in the surrounding area of D1 and D2.The accelerations occurred on D1 and D2 after January 2017 might correspond to land conversions in Fig.12(c) and (d).

Figure 12 :
Figure 12: (a) Subsidence rate of Dongxihu area during 2015 -2019.The rectangle represent the location of (c) and (d).(b) is the 235
). Signs of construction were observed from the optical image acquired at February 2016 in Fig. S2(b).Acceleration from the time series subsidence was also observed in Fig. 11.Similar behaviours can be observed also from QS2, which is about 500 m away from the Shisan-Jiefang and Shisi-Jiefang.Interpretation from the optical image acquired in May 2017 in Fig. S2(c), the housing demolition is almost finished.Constructions started since May 2018 and buildings are identified at Shisan-Jiefang from the optical image acquired in October 2018 Fig. S2(d).Accelerations can be clearly observed by InSAR measurements in Fig. 11 during the deep foundation pit dewatering.Meanwhile, the first terrace and second terrace EGSs are composed of typical binary structural stratums (Table

Figure 14 :
Figure 14: (a) Nonlinear subsidence of HH1 and QL1, (b) Water level of Yangtze River and rainfall.