Surface Deformation Evolution in the Pearl River Delta Between 2006 and 2011 Estimated from the ALOS1/PALSAR

This study monitors the land subsidence of the whole Pearl River Delta (PRD) (area: ~40,000 km 2 ) in China using the ALOS1/PALSAR data (2006-2011) and the SBAS-InSAR method. We also analyze the relationship between the subsidence and the coastline change, river distribution, geological survey data as well as the geomorphology data. The results show that (1) the land subsidence with the average velocity of 50 mm/year occurred in the low elevation area in the front part of the delta and the coastal area; (2) the areas along rivers also suffered from surface subsidence, due to the fragile geological structure and frequent human interference; (3) the geological evolution is the intrinsic factor of the surface subsidence in the PRD, but human interference (reclamation, ground water extraction and urban construction) extend the subsiding area and accelerate the subsiding rate. This study would provide technique supports for the surface deformation and disaster warning in delta regions.


Introduction
Deltas, where sea and land interact, have abundant natural resources, superior natural environment and geographical locations, so they accommodate the most intensive human activities. However, most deltas suffer from geological disasters caused by land subsidence. Land subsidence damages urban infrastructure, such as buildings, railways, highways, subways and underground facilities (water supply pipes, gas and electrical installations), hindering the sustainable development in economy and environment (Ericson et  The interferometric synthetic aperture radar (InSAR) technique has been widely used for monitoring surface subsidence in delta areas, because of its wide spatial coverage and continuous temporal coverage (Higgins et al. 2014). The ENVISAT/ASAR, ALOS1/PALSAR ,Sentinel-1 and COSMO-SkyMed data have been used to monitor the surface deformation of the whole PRD or parts of the PRD, and the results suggested that most local land subsidence was mainly caused by human activities (ground water extraction, engineering construction et al.) (Wang et Du et al. 2020). These studies all attribute the land subsidence to human activities, which may be not inclusive. Actually, natural evolution factors (such as river distribution, coastline change and geomorphology) may also lead to the land subsidence. However, the relationship between those factors and the surface subsidence are still unclear.
In this study, we apply the small baseline subset InSAR (SBAS-InSAR) method to map the surface deformation over the entire PRD (area: ~40,000 km 2 ) using 162 ALOS1/PALSAR images from 3 paths and 8 frames acquired between 2006 and 2011. The deformation results are veri ed by cross-validation.
Then the temporal-spatial distribution characteristics of the deformation in the PRD are analyzed. Finally, we analyze the relationship between the subsidence and natural evolution factors as well as human interferences, and discuss the surface deformation pattern of the PRD and the implications.

Geological background
The PRD is located in the south-central part of Guangdong province, bordering the South China Sea ( Fig. 1). It is a composite delta plain formed by the alluvium from the Xijiang River, Beijiang River and Dongjiang River (Huang and Zhang 2004). It covers an area of about 42, 200 km 2 from 112.5° to 114.25°E and 21.83° to 23.6°N. This delta is high in the north and low in the south, with small undulations. In uenced by the subtropical monsoon climate, the PRD has abundant rainfall and is often struck by typhoons and rainstorms. More than 300 river sections, with a total length of about 1740 km and an area of over 9,750 km 2 , form a dense and complex network (Huang and Zhang 2004).
Quaternary loose soil (Q4), bedrock (β), hidden karst (C) and soft soil form the geological structure of this area (Fig. 2). The bedrock is mainly composed of granite and partially of metamorphic sandstone, clastic rocks and basalt. The concealed karst, composed of limestone and shale, concentrates in the Guanghua Basin. The Quaternary loose soil and soft soil form alluvial layers, which are thick in the front part (such as Zhongshan city and Nansha district of Guangzhou city) and thin in the rear part of the delta. The soft soil layer was formed late, mostly from the Pleistocene to the Holocene (Zong et al. 2009). The thickness of most parts of the layer is between 5 and 36 m, but the maximum is 55 m. The soft soil layer mostly distributes along the river and in the at terrain of the delta (Li et al. 1991).

River and coastline evolution of the PRD
The rich rainfall and abundant river ow carry sediments from the upper reaches to the lower reaches. The PRD is at, so 20%~50% of the total sediment transport deposits in the river channel (Hu et al. 2011).
Since the Tang dynasty, the population of the PRD started soaring. The agricultural development brought by the population growth accelerated the change of land use and led to serious soil and water loss. The

SBAS-InSAR data processing
For a differential InSAR pair acquired at moments t A and t B , the unwrapped phase can be written as where δφ obs is the phase of one pixel in differential interferogram, Δd is the line of sight (LOS) deformation of a pixel between two images, Δz is the topographic error, Δφ atm is the atmospheric delay phase, Δ n j is the phase noise. Eq. (1) can be written in a matrix form as follows: where v is the deformation rate. Following the SBAS-InSAR procedures (Berardino et al. 2002), we solve Eq. (2) by the least squares method, and obtain the linear deformation rate and the topographic error of each pixel. Subtracting topographic error and linear deformation phase from each differential interferogram, we get the residual phase. The main components of the residual phase are nonlinear deformation, atmospheric delay and noise. Taking the residual phase as the observation value, we obtain the time series phase of each SAR image by the singular value decomposition (SVD) method. According to the characteristics of nonlinear deformation, atmospheric delay and noise in the space-time domain, we use the temporal high pass and spatial low pass ltering to obtain the atmospheric phase of each SAR image. Then, atmospheric errors are subtracted from each interferogram. Finally, the SVD method is used to calculate time series of each pixel. Figure 3 shows the owchart of the SBAS method.
In this study, the L-band ALOS1/PALSAR data from 3 tracks and 8 frames are used to obtain the deformation in the PRD occurred between 2006 and 2011 (see Fig. 1, and Appendix A). The data is processed with GAMMA software. The multi-look operation of 6 × 16 pixels in the range and azimuth directions is adopted to reduce phase noise and calculation before phase unwrapping. The image resolution after multi-look is about 50 m × 50 m. The interferometric pairs with the temporal baseline less than 700 days and spatial baseline less than 1000 m are selected to reduce the decorrelation. The topographic phase is removed using the 1-arcsecond (~ 30 m) Shuttle Radar Topography Mission DEM (Jarvis et al. 2008). A modi ed Goldstein lter is used to reduce noise phases (Li et al. 2008) and a mask with a coherence less than 0.4 is adopted to weaken the effect of low coherence regions (sea and isolated islands, etc.). Finally, the phase unwrapping is performed by the minimum cost ow algorithm (Mario 1998). The polynomial tting is used to remove the long wave-length orbital errors in the unwrapped phase.

Surface subsidence results
Using the method described in Sect. 3.1, we obtain the deformation rate map of the PRD between 2006 and 2011 ( Fig. 4). Figure 4a presents the average deformation rate obtained by the ALOS1/PALSAR image. The hinterland is stable and the subsidence is mainly distributed in the front part of the delta and the coastal area. The average subsidence rate ranges from 10 to 40 mm/y, with the maximum deformation rate of > 70 mm/y. Some regions have sporadically distributed subsidence with a rate between 10 to 30 mm/y, such as the northeast part of Guangzhou city and the southeastern part of Zhaoqing city (outlined by the red dotted line in Fig. 4).
Two large subsidence cones (S-1 and S-2) are observed in the front part of the delta (see Fig. 4b). The most serious land subsidence is in the junction area of Zhongshan city and Jiangmen city (such as Jianghai district of Jiangmen city and Henglan town of Zhongshan city), with an area of 160 km 2 . The subsidence rate ranges from 30 to 50 mm/y, even > 70 mm/y in some parts. The average cumulative deformation of the region is about 35 cm from 2006 to 2011 (Fig. 4d). S-2 covers Jinwan district of Zhuhai city and Tanzhou town of Zhongshan city, with an area of 230 km 2 and a subsidence rate of 30-50 mm/y. Between the two large subsidence cones, there are some areas with smaller subsidence rates, such as Doumen district of Zhuhai city and Xinhui district of Jiangmen city.
No large subsidence cone is found in the coastal areas. The subsidence areas scatter on the two sides of the Lingding Bay (see Fig. 4c), such as Nanlang district of Guangzhou city and Nanlang town of Zhongshan city and Tangjia town of Zhuhai city. The subsidence rate is between 20 and 40 mm/y, and is over 50 mm/y in some regions. The subsidence rates in the coastal areas of Dongguan and Shenzhen are smaller. For example, the subsidence rates of Chang'an town in Dongguan and Bao'an district in Shenzhen are 10 and 25 mm/y, respectively.

Results veri cation
In order to verify our results, we compared the observations of Area1 and Area2 in Path 461 and Path 460, which are the overlapping areas of those two paths at the same time period (Fig. 4(a)). Area1 is stable but Area2 has obvious subsidence, so they can verify the results of different subsidence levels.
As Fig. 5 shows, the velocity difference in Area1 is between − 8 and 5 mm/y, and the overall mean value is -1.3 mm/y with a variance of 2.2 mm/y. In Area2, the velocity difference is between − 10 ~ 6 mm/y, the mean value is -0.4 mm/y, and the variance is 3.8 mm/y. Because of decorrelation, atmospheric delays and different incidence angles, the results of the overlapping areas of adjacent paths, like Area1 and Area2, are slightly different, so the uncertainty (2.2 mm/yr and 3.8 mm/yr) between different paths is reasonable. To further verify the results, we did eld investigations in some cities of the PRD, including Jiangmen city, Zhuhai city and Guangzhou city (see Fig. 6). The led observation is perfectly consistent with the model results. Superimposing the river distribution map of the PRD on the deformation rate map (see Fig. 7a), we also nd that the subsidence areas are distributed along rivers, especially the front part of the PRD, where the river networks are dense (Chen et al. 2012). Moreover, the map of the coastline changes over the past 6,000 years (see Fig. 7b) show that almost all the subsidence areas are distributed in the regions formed in the past 1000 years. The coastline of the PRD continues pushing outward as silt keeps on depositing at the junction of land and sea. The geological structure formed by subsidence, however, is inherently unstable. The PRD is mainly composed of soft soil layers with high water content and compressibility, poor permeability, low bearing capacity and shear strength (Gebrechael et al. 2018), which easily leads to surface deformation. Most of the soft soil layers formed by sediment deposition are distributed along rivers, so subsidence usually occurs along rivers and in the front part of the delta.
The formation time of the coastline also has in uences on the stability of the ground structure, as newly formed soil has a shorter consolidation period and may bring larger subsidence. The land subsidence of the PRD has a high correlation with the distribution Quaternary sediments (Chen et al. 2012). This pattern has been validated in various delta regions of the world without man interventions.

Geological structure
Some areas in the front part of the PRD has no signi cant settlement (see Fig. 7). These areas are the bedrock distribution areas (β) (see Fig. 8a). Before the modern delta were formed (6,000 years ago), these bedrock distribution areas were isolated islands, and most of them were composed of bare rocks. As the coastline moved outward, those islands gradually merged into the land. The geological structures of these regions are stable. As Fig. 8a and 7b show, there is a strong positive correlation between the thickness of the soft soil and the deformation rate. Similar results were got in the research on the Ganges-Brahmaputra Delta (Higgins et al. 2014). The relationship between the compressible soil layer thickness and subsidence rate in the Tiber River delta was also investigated with the borehole data and the groundshaped variables obtained from the PS-InSAR technique (del Ventisette et al. 2015). Therefore, the soft soil thickness can be inverted by the correlation between surface deformation and drilling data.
Soft soil is mainly distributed along rivers and low-lying areas of the delta, where the erosion caused by sea level rise is serious. In Fig. 8c, the green area is lower than the average sea level, and it accounts for more than 1/3 of the entire onshore delta. The subsidence in the area below the sea level is outlined by the yellow line in Fig. 8d. The comparison between Fig. 8c and Fig. 8d shows that the land subsidence in the front and middle parts of the PRD concentrates in the area below the mean sea level. The annual average tidal range of the Pearl River Estuary is 1.0-1.7 m, and the load generated by tides and storm surges can cause settlements of a few millimeters to centimeters. The low-lying areas, affected by tides, rain, sea level rise, storms and oods, are also the areas accommodating severe subsidence disasters.
According to the Intergovernmental Panel on Climate Change (IPCC), the sea level had risen 0.19 m (1.7 mm/y) between 1901 and 2010. And the rise is accelerating (Zhang and Church 2012). These lowlying areas would suffer from both seawater erosion and surface subsidence in the future.

Human activities 4.2.1 Reclamation
In the 1960s, the large-scale reclamation in the PRD rapidly changed the shape of the coastline. The reclamation rate between 1970 and 1995 reached 21 km 2 /y and then slowed down in the past 10 years. We collect the coastline data from 1850 to 2015 with a time span of 165 years, and depict the changes in Fig. 9. The blue line is the PRD coastline in 1850, and the red lines are the coastlines of different periods.
In the past 165 years, the coastline and landforms have undergone tremendous changes. A large area of land has been reclaimed from waters and tidal ats. In particular, about 1,258 km 2 water area of the estuary bay has been reclaimed in the past 45 years (from 1970 to 2015). As the shoreline gradually moved towards the sea, more than 20 islands gradually merged into the mainland. By 2015, almost all the tidal ats in 1850 had been reclaimed into land.
As shown in Fig. 9, the reclaimed areas are also the most severely settlement area (subsidence rate: 50 mm/y), including the coastal areas of Zhuhai, Guangzhou, Dongguan and Shenzhen. Most of these areas were reclaimed in the last 30 years (Xu et al. 2016; Huang and Zhang 2004), so the geological structure is unstable. Furthermore, the land ll materials are sediments, which causes subsidence by selfconsolidation.
We select two areas, Nansha district of Guangzhou city and Jinwan district of Zhuhai city (outlined by the yellow dotted line in Fig. 9) to analyze the relationship between the subsidence and reclamation in detail. As shown in Fig. 10, these two areas were reclaimed from the sea a few decades ago. In 1984 (see Fig. 10a), Longxue and Hengmen were small islands, but they gradually expanded outwards and formed the basic shape of the current island in 2016 (see Fig. 10b-e). The area of Longxue island cofferdam had increased about 150 times in 35 years. The newly reclaimed regions have serious subsidence, and the maximum subsidence rate is 75 mm/y (see Fig. 10f). Figure 10g-k show the variation of the reclaimed region in Jinwan district of Zhuhai city, where two small subsidence cones are obvious (see Fig. 10l).
Gaolan island, which was an island before 2004, became connected to the land in 2010 (see Fig. 10j). The connection area accommodates the subsidence cones, where the maximum subsidence is 87 mm/y (see Fig. 10l).
Many countries, such as the Netherlands (Hoeksema 2007), the United Kingdom (Wali 1975), Japan (Suzuki 2003) and Singapore (Douglas and Lawson 2003;Xu et al. 2016), have reclaimed land from the sea for the coastal industries and agricultural development, as well as for storm defense. Surface deformation monitoring in the Nile Delta and the Tiber Delta indicates that the reclaimed area is the most severely subsiding area (Stramondo et al. 2008;Marriner et al. 2012). So the long-term monitoring of the reclaimed area in the PRD is necessary to ensure the safety of the construction there.

Other human interventions
The PRD region is an economic center of south China, where densely distributed factories, buildings, farmlands, and sheries. The natural evolution and reclamation are the main factors for land subsidence in the PRD, but the land subsidence rate in some areas is as high as 70 mm/y, which can not only be caused by natural factors. We took the subsidence cone S-1(see Fig. 4b) in Jianghai District of Jiangmen City as an example to analyze the relationship between human activities and surface subsidence. As shown in Fig. 11b, the maximum subsidence rate is 80 mm/y.
The eld investigation found that sh ponds are widely distributed in the coastal area of Jianghai (Fig. 11a), where large area land subsidence occurred. Aquaculture needs great amount of freshwater. Aggressive extraction of groundwater leads to the rapid drop of the groundwater level and causes subsidence, even settlement funnels. Settlement funnel S-1 is 160 km 2 and is sinking 50 mm/y. Two pro les along the settlement funnel S-1 (Fig. 11b), ab and cd, are used to analyze the settlement characteristics and their average settlement rates are shown in Fig. 12.
The subsidence caused by excessive groundwater extraction is common in delta areas. In the Mekong Delta, the average land subsidence rate has reached 16 mm/y as more than 250,000 hectares of land has been converted from rice elds to shrimp ponds since 2000 (Erban et al. 2014). Oil and groundwater extraction play a key role in the land subsidence in the modern Yellow River Delta, and the land subsidence caused by aquaculture is up to 250 mm/y (Higgins et al. 2013). The land subsidence caused by human activities usually has a small involved area, about a few square kilometers to tens of square kilometers. The area of the subsidence funnel S-1 in the PRD is up to 160 km 2 , such magnitude is of rare occurrence in other deltas.
Natural factors coupled with human activities cause the subsidence in the PRD. The evolution of the PRD determines the location of the subsidence areas, and human activities have become the direct cause of the settlement. Therefore, in addition to the subsidence caused by natural factors, we should pay more attention to the subsidence caused by human activities.

Conclusions
We map the surface deformation of the entire PRD using the ALOS1/PALSAR images from 2006 to 2011.
The results show that the areas with serious subsidence concentrate in the front part of the delta and the coastal areas. The two large settlement funnels in the front part of the PRD have the average settlement rates of 30-50 mm/y, and the maximum settlement rate of larger than 70 mm/y. The surface deformation map generated from the InSAR observation shows that: (1) the subsidence area is obviously distributed along rivers; (2) the natural evolution would be the intrinsic factor for the land subsidence of the PRD; (3) the land subsidence in the PRD is signi cantly related to human activities (such as land reclamation, pumping groundwater for aquaculture). Due to the natural and human factors, the settlement in the PRD region will continue. These ndings have guiding contributions to disaster prevention and social construction in the PRD.      The histogram of the difference of the average deformation rate of (a) Area-1 and (b) Area-2 derived from different frames.

Figure 6
Photos of the led investigation in the areas denoted by the yellow stars in Figure 4. (a) Xinhui port in Xinhui district, Jiangmen city (P-1 in Figure 4). The ground subsidence is about 5 cm under the steps; (b) the Kangdelai Medical city in Jinwan district, Zhuhai city (located at P-2 in Figure 4). The steps have separated from the building; (c) a building in Honghai village, Wanhaisha town, Nansha district, Guangzhou (P-3 in Figure 4). The underground pipeline has exposed because of the ground subsidence; (d) a house in Sanban village, Hongqi town, Jinwan district, Zhuhai city (P-4 in Figure 4). The windows have submerged into the ground.