Evaluating Potential Ground Subsidence Geo-Hazard of Xiamen Xiang’an New Airport on Reclaimed Land by SAR Interferometry

: The land reclaimed from the seaside may have a long-term subsidence trend, which poses a potential geohazard in the future land use. Xiamen Xiang’an New Airport (XXNA) is built on reclaimed land since 2016. Based on the spaceborne Sentinel-1 data between January 2018 to April 2019 and the time series interferometric synthetic aperture radar (InSAR) technique, this paper analyzed the reclaimed land subsidence evolution at XXNA in this period. InSAR measurements show that XXNA is su ﬀ ering from severe subsidence, mainly in three regions because of the earth and sand compacting. By analyzing the spatial subsidence characterizations of the main subsiding areas combined with historical land reclamation and future land use planning, we ﬁnd the potential threat of subsidence to future land use. Correlation between subsidence and the period of reclamation was found, indicating that the consolidation and compression in dredger ﬁll is the main cause of subsidence. By combining subsidence monitoring results with di ﬀ erent land use types and adopting the Expectation (Ex) and Entropy (En) methods, we analyzed the key area with potential subsidence geo-hazard. This work shows that with SAR interferometry, it is possible to ﬁnd the large area ground subsidence in the airport reclaimed area. The areas with potential subsidence geo-hazards are consistent with the deep reclaimed earth, which means high subsidence risk in the future.


Introduction
With the urbanization and economic development in coastal areas, land resources are becoming increasingly scarce. To meet the development needs of coastal cities, humans have to reclaim land from the seaside all over the world, such as Shenzhen (China), Shanghai (China), Urayasu (Japan), New Orleans (USA), and Venice (Italy). All of those cities [1][2][3][4][5] have suffered from the land subsidence hazards due to long-time and gradual soil compaction on reclaimed land, which would cause damage to important infrastructures (e.g., airports, roads, buildings, and metros) and pose a great threat to the environment and public safety [6]. Therefore, monitoring, analyzing, and controlling the subsidence Sustainability 2020, 12, 6991 2 of 17 on reclaimed land is a crucial task not only to ensure the safety of citizens' life and property, but also to support the sustainable development of coastal areas.
Xiamen Xiang'an New Airport (XXNA), one of the most important international airports in China, is built on Dadeng Island in Xiamen, Fujian Province. XXNA was scheduled to start construction in 2016, enter service in 2020, and complete construction by 2045 [7]. The planning area of Dadeng Island reaches 43 km 2 , and about 65% of the land is reclaimed from the sea. According to previous studies [7,8], due to the soil consolidation of the underlying unconsolidated marine sediments [9], XXNA is undergoing significant and non-uniform subsidence on reclaimed land during its construction. Therefore, it is important and necessary to continuously monitor ground subsidence at XXNA and determine potential subsidence damage during its construction and for future airport operation.
Traditional in-situ measurements, such as leveling and Global Navigation Satellite System (GNSS), have been widely employed in ground subsidence monitoring with high accuracy. However, they cannot acquire continuous and dense ground subsidence measurements with a large-scale coverage at a low cost [10,11]. Over the past decades, interferometric synthetic aperture radar (InSAR) technique has been proven to be a powerful geodetic technique to measure surface displacements with high-accuracy. It has been successfully applied in city subsidence detection [12][13][14], landslides monitoring [15][16][17], earthquake displacements measurements [18][19][20], etc. With the advantage of weather independence, wide coverage, high resolution, and high accuracy [21][22][23], InSAR technique has made many achievements in airport and reclaimed land subsidence monitoring. To monitor the subsidence in airports, Zhao et al. [24] and Ding et al. [25] investigated residual settlement of the reclaimed foundation at Hong Kong International Airport by InSAR technique. Jiang et al. [26] employed high-resolution TerraSAR-X images to study the subsidence characteristics and driving forces of subsidence on Pudong International Airport from 2011 to 2012. Marshall et al. [27] applied a multi-temporal InSAR method to investigate the subsidence on Kuala Lumpur International Airport with Sentinel-1 images from 2014 to 2017. Dai et al. [28] used high-resolution Cosmo-SkyMed SAR data from 2017 to 2019 at Beijing Capital International Airport to diagnose the subsidence hazard. In terms of monitoring reclaimed land subsidence, Xu et al. [1] applied multi-temporal InSAR method to investigate coastal reclaimed land subsidence in Shenzhen by an ENVISAT ASAR ascending and descending orbits dataset both acquired from 2007 to 2010. Yang et al. [2] determined coastal subsidence in Shanghai combined with historical reclamation activities by TerraSAR-X, ALOS PALSAR, and ENVISAT ASAR images from 2007 to 2010. Aimaiti et al. [3] investigated land reclamation in the coastal city of Urayasu by combining the ERS-1/-2 and ALOS-2 PALSAR-2 datasets from 1993 to 2017. Liu et al. [8] employed SBAS-InSAR technique to determine the temporal and spatial evolution of ground settlement at XXNA combined with land reclamation (with Sentinel-1 images from 2015 to 2018). Although previous studies have shown good examples on airport or reclaimed land subsidence monitoring by InSAR technique, the relationship between subsidence on reclaimed land and future land use has not yet been qualitatively analyzed. The subsidence risk of XXNA would have great potential risk for the safe operation of the whole airport, which needs an advanced and continuous high-accuracy subsidence monitoring.
In this paper, the time series InSAR technique was utilized to monitor reclaimed land subsidence at XXNA based on 36 ascending Sentinel-1 SAR images acquired from January 2018 to April 2019. The main subsiding areas were found, and spatiotemporal subsidence characteristics were analyzed by combining with historical land reclamation and land use planning. Finally, we discussed potential risks of subsidence at XXNA by different quantitative statistical indicators.

Location of XXNA and Land Use Planning
XXNA is located at the southeast of the coastal region in China. It belongs to Xiang'an district, Xiamen, China, 25 km to the west of the heart of Xiamen (Figure 1). XXNA will be built in the east of Dadeng Island (Figure 1), aiming to be as regional hub airport and major gateway airport. In the short-term target plan (target year 2025), the construction of two long-distance parallel runways is planned, and the annual passenger throughput is expected to reach 45 million passengers, and cargo and mail throughput of 750,000 tons. The long-term goal (target 2045) is to add a shortrange parallel runway to each of the two long-distance parallel runways. By then, XXNA will have four runways, which can meet the annual passenger throughput of 70 million passengers, and cargo and mail throughput of 1.6 million tons of transportation requirements.   According to the land use planning of XXNA in Figure 2 (adapted from Zhao [7]), the total planned area on Dadeng Island is about 43 km 2 , including eighteen different land use types. The area for airport is about 16.30 km 2 , which accounts for 37% of all planning area on Dadeng Island. The construction plan for XXNA is divided into two phases: short-term goals and long-term goals [7,29]. In the short-term target plan (target year 2025), the construction of two long-distance parallel runways is planned, and the annual passenger throughput is expected to reach 45 million passengers, and cargo and mail throughput of 750,000 tons. The long-term goal (target 2045) is to add a short-range parallel runway to each of the two long-distance parallel runways. By then, XXNA will have four runways, which can meet the annual passenger throughput of 70 million passengers, and cargo and mail throughput of 1.6 million tons of transportation requirements. According to the land use planning of XXNA in Figure 2 (adapted from Zhao [7]), the total planned area on Dadeng Island is about 43 km², including eighteen different land use types. The area for airport is about 16.30 km 2 , which accounts for 37% of all planning area on Dadeng Island. The construction plan for XXNA is divided into two phases: short-term goals and long-term goals [7,29]. In the short-term target plan (target year 2025), the construction of two long-distance parallel runways is planned, and the annual passenger throughput is expected to reach 45 million passengers, and cargo and mail throughput of 750,000 tons. The long-term goal (target 2045) is to add a shortrange parallel runway to each of the two long-distance parallel runways. By then, XXNA will have four runways, which can meet the annual passenger throughput of 70 million passengers, and cargo and mail throughput of 1.6 million tons of transportation requirements.

Historical Land Reclamation
XXNA on Dadeng Island was reclaimed from the sea to meet the needs of urban development. Remote sensing methods provide us with an efficient way for land change detection [30,31]. To investigate the historical land reclamation on Dadeng Island, Landsat-8 satellite images and Google Earth historical images from 2013 to 2019 were utilized to assist our analysis. The results of the temporal and spatial evolution of land reclamation are shown in Figure 3, 16 km 2 of land, which accounts for 36% of all, has been reclaimed from the sea from 2013 to 2019. When the reclamation project is complete, reclaimed land will reach 29 km 2 (accounts for 65% of all). As shown in Figure 3a, the whole reclamation work is mainly in the southeast of Dadeng Island, the northwest and southwest also have wide reclaimed area. Additionally, according to previous studies [7,8], the reclaimed land at XXNA is mainly composed of sand-mixed silt (Q 4 m ), silt and silt mixed with sand (Q 4 m ), silty clay and clay (Q 3 al ), silty sand (Q 3 al ), and residual clay (Q el ) (The Q 4 m , Q 3 al , and Q el indicate the quaternary marine sediments, the quaternary alluvium sediments, and the quaternary eluvial sediment, respectively). When there is a heap load on the surface, it is prone to obvious subsidence or non-uniform settlement on those land types.

Historical Land Reclamation
XXNA on Dadeng Island was reclaimed from the sea to meet the needs of urban development. Remote sensing methods provide us with an efficient way for land change detection [30,31]. To investigate the historical land reclamation on Dadeng Island, Landsat-8 satellite images and Google Earth historical images from 2013 to 2019 were utilized to assist our analysis. The results of the temporal and spatial evolution of land reclamation are shown in Figure 3, 16 km 2 of land, which accounts for 36% of all, has been reclaimed from the sea from 2013 to 2019. When the reclamation project is complete, reclaimed land will reach 29 km 2 (accounts for 65% of all). As shown in Figure  3a, the whole reclamation work is mainly in the southeast of Dadeng Island, the northwest and southwest also have wide reclaimed area. Additionally, according to previous studies [7,8], the reclaimed land at XXNA is mainly composed of sand-mixed silt (Q4 m ), silt and silt mixed with sand (Q4 m ), silty clay and clay (Q3 al ), silty sand (Q3 al ), and residual clay (Q el ) (The Q4 m , Q3 al , and Q el indicate the quaternary marine sediments, the quaternary alluvium sediments, and the quaternary eluvial sediment, respectively). When there is a heap load on the surface, it is prone to obvious subsidence or non-uniform settlement on those land types.

Datasets
Based on the previous monitoring [8], as the reclaimed land surface is still undergoing substantial change, long time span images will cause SAR image decorrelation, causing the newly reclaimed area to be impossible to monitor. Therefore, in this study, we employed 36 ascending Sentinel-1A TOPS mode Interferometric Wide Swath (IWS) images, acquired from 3 January 2018 to 16 April 2019, to ensure high correlation and estimate land reclamation evolution and time series subsidence. The Sentinel-1 satellites developed by the Europe Space Agency (ESA) of Copernicus

Datasets
Based on the previous monitoring [8], as the reclaimed land surface is still undergoing substantial change, long time span images will cause SAR image decorrelation, causing the newly reclaimed area to be impossible to monitor. Therefore, in this study, we employed 36 ascending Sentinel-1A TOPS mode Interferometric Wide Swath (IWS) images, acquired from 3 January 2018 to 16 April 2019, to ensure high correlation and estimate land reclamation evolution and time series subsidence. The Sentinel-1 satellites developed by the Europe Space Agency (ESA) of Copernicus Program conduct C-band SAR imaging with a wavelength of 5.6 cm. Compared with other SAR missions, Sentinel-1 has the characteristics of short revisit time (twelve days for one satellite, and six days for two satellites constellation) and short baselines (i.e., orbital separations). Detailed information and main parameters regarding to the Sentinel-1 data are shown in Table 1. We set the image 20180328 as master image for co-registration with other images. With a threshold of temporal baseline below 60 days and spatial perpendicular baseline below 105 m, a total of 137 interferograms were generated. Detailed spatial and temporal baselines of these Sentinel-1 datasets are shown in Figure 4. ALOS World 3D 30 m Digital Elevation Model (DEM) was used to remove the topography-related phase and assist the geocode processing. Program conduct C-band SAR imaging with a wavelength of 5.6 cm. Compared with other SAR missions, Sentinel-1 has the characteristics of short revisit time (twelve days for one satellite, and six days for two satellites constellation) and short baselines (i.e., orbital separations). Detailed information and main parameters regarding to the Sentinel-1 data are shown in Table 1.
We set the image 20180328 as master image for co-registration with other images. With a threshold of temporal baseline below 60 days and spatial perpendicular baseline below 105 m, a total of 137 interferograms were generated. Detailed spatial and temporal baselines of these Sentinel-1 datasets are shown in Figure 4. ALOS World 3D 30 m Digital Elevation Model (DEM) was used to remove the topography-related phase and assist the geocode processing.

SBAS-InSAR Method
The SARscape Modules for ENVI software suite (HARRIS Geospatial Solutions, Broomfield, CO, USA) were used to perform the Small Baseline Subset InSAR (SBAS-InSAR) processing [32,33]. The main steps were as follows. There are N+1 SLC (Single Look Complex) images acquired according to

SBAS-InSAR Method
The SARscape Modules for ENVI software suite (HARRIS Geospatial Solutions, Broomfield, CO, USA) were used to perform the Small Baseline Subset InSAR (SBAS-InSAR) processing [32,33]. The main steps were as follows. There are N+1 SLC (Single Look Complex) images acquired according to time series t 0 . . . t i . . . t N . The first step of the SBAS method was "connection graph", which selected SAR image pairs with spatiotemporal baseline constraints. M interferograms were generated after the removal of the topographic phase based on orbital information and external DEM data in the step of "differential interferogram generation". Multi-looking in range and azimuth direction and the Goldstein filtering method was used to improve the signal/noise ratio of the interferograms [34]. Then followed the completion of phase unwrapping of the wrapped phase signal calculated from the selected interferograms by Minimum Cost Flow (MCF) method [35][36][37]. Several interferogram pairs with unwrapped phase error or low coherence were be discarded [38]. The next step was "refinement and re-flattening". By selecting and refining the stable Ground Control Point (GCP), the residual phase content and phase ramps were calculated to correct the unwrapped phase [39]. Then came the key step of "first estimation of residual height and the displacement velocity". In this step, assuming that the time t 0 is used as the reference time, and the subsidence of the study area is regarded as zero, the phase of a pixel in the j-th interference image relative to the unwrapping start point can be expressed as follows: where, 1 ≤ j ≤ M, δϕ j def is subsidence phase, δϕ j topo represents topography phase, δϕ j atm represents atmosphere phase, and δϕ j noise represents noise phase. Those could be further expressed as follows: where, λ is the radar wavelength, θ represents the radar angle, ∆h accounts for the error of DEM, r accounts for the slant range. d(t 1 ) and d(t 2 ) represent the cumulative amount of subsidence compared to the initial time (t 0 ) in the direction of the radar sight. We generally assume that d(t 0 ) = 0.
In the SBAS method, it is assumed that the surface subsidence in the interval between the acquisition of two adjacent SAR images conforms to the linear cumulative trend [39], and the solution of the phase time series is converted into the solution of the phase change rate. It could be expressed as follows: Then δϕ j def could be expressed as: In Equation (6), matrix B is a coefficient matrix, and each row corresponds to an interference pair. Among the elements of the matrix, the main image coefficient is 1, the secondary image coefficient is −1, and the other image coefficients are 0.
Therefore, without considering atmospheric phase and noise, the joint expressions (1)-(4) can be shown as: If the interference network forms a unique subset and is fully connected, the least square method can be used to solve the subsidence rate and the elevation correction number based on equation (5). When the matrix has a rank deficit, it can be solved by the Singular Value Decomposition (SVD) method. After the first estimation is completed, the atmospheric phase is estimated by atmospheric filtering. Then, the step of "second estimation of residual height and the displacement velocity" will remove the atmospheric phase and estimate as first estimation. Finally, the rates of land subsidence, time series of subsidence, and residual height are derived and mapped across the study area. To evaluate the average precision of the displacement velocity, the standard deviation value of the measurement was derived based on the parameters' coherence and wavelength. The formula used for the precision calculation is: where γ is the interferometric coherence. The higher this value, the lower the measurement precision. The flowchart of SBAS-InSAR processing was presented in Figure 5.
Sustainability 2020, 17, x FOR PEER REVIEW 7 of 17 method. After the first estimation is completed, the atmospheric phase is estimated by atmospheric filtering. Then, the step of "second estimation of residual height and the displacement velocity" will remove the atmospheric phase and estimate as first estimation. Finally, the rates of land subsidence, time series of subsidence, and residual height are derived and mapped across the study area. To evaluate the average precision of the displacement velocity, the standard deviation value of the measurement was derived based on the parameters' coherence and wavelength. The formula used for the precision calculation is: Where γ is the interferometric coherence.
The higher this value, the lower the measurement precision. The flowchart of SBAS-InSAR processing was presented in Figure 5.

Expectation and Entropy Methods
The Expectation (Ex) and Entropy (En) methods were proposed in the cloud model, which represents the center of gravity of the cloud and the measure of the fuzziness of the concept, respectively [40,41]. For land subsidence research, Ex and En can be used to represent the overall level and the degree of non-uniformity land subsidence, respectively [42,43]. The high value of Ex means more significant land subsidence, and the high value of En represents more serious inhomogeneity of land subsidence.
The values of Ex and En can be calculated as follow: Where xi is the value of the high-coherence points and N is the number of high-coherence points achieved by InSAR method in study area. Then En can be calculated by: In previous studies [42,43], the Expectation and Entropy methods have been used to characterize the spatiotemporal difference of the displacement along subway and railway. In this paper, the Expectation and Entropy methods were used to analyze the overall level and the degree of nonuniform land subsidence in different land use types.

Expectation and Entropy Methods
The Expectation (Ex) and Entropy (En) methods were proposed in the cloud model, which represents the center of gravity of the cloud and the measure of the fuzziness of the concept, respectively [40,41]. For land subsidence research, Ex and En can be used to represent the overall level and the degree of non-uniformity land subsidence, respectively [42,43]. The high value of Ex means more significant land subsidence, and the high value of En represents more serious inhomogeneity of land subsidence.
The values of Ex and En can be calculated as follow: where x i is the value of the high-coherence points and N is the number of high-coherence points achieved by InSAR method in study area. Then En can be calculated by: In previous studies [42,43], the Expectation and Entropy methods have been used to characterize the spatiotemporal difference of the displacement along subway and railway. In this paper, the Expectation and Entropy methods were used to analyze the overall level and the degree of non-uniform land subsidence in different land use types.

Mean Velocity Subsidence and the Time Series Subsidence
With Sentinel-1 ascending datasets from January 2018 to April 2019, the mean velocity subsidence map on Dadeng Island, as shown in Figure 6a, was generated by SBAS-InSAR method. The negative values mean that the ground is moving away from the satellite. According to previous studies [1,8], the vertical displacement dominates ground deformation for subsidence on reclaimed land. When the horizontal displacement is much lower than the vertical displacement, the vertical subsidence dominates the LOS displacements [28]. As XXNA was in flat topography, the displacements acquired from the InSAR method were reasonably interpreted as vertical subsidence in this study. As we can see in Figure 6a, severe subsidence can be observed at XXNA (rectangle A in Figure 6a) and the subsidence rates can reach up to 192 mm/year, which is the position with the largest subsidence rate in the entire study area. In Qiaotou (rectangle B in Figure 6a) and Dengqi (rectangle C in Figure 6a), the maximum subsidence rates can reach 169 mm/year and 156 mm/year, respectively. Except for these three regions mentioned, the rest of Dadeng Island is relatively stable. Figure 6b presents the error map of mean velocity results calculated by the standard deviation value as shown in equation 7. It can be seen that there is high precision in most parts of study area (below ±5 mm/year), and the precision in a very small part of study area (with red color in Figure 6b) is estimated as 5~12 mm/year.

Mean Velocity Subsidence and the Time Series Subsidence
With Sentinel-1 ascending datasets from January 2018 to April 2019, the mean velocity subsidence map on Dadeng Island, as shown in Figure 6a, was generated by SBAS-InSAR method.
The negative values mean that the ground is moving away from the satellite. According to previous studies [1,8], the vertical displacement dominates ground deformation for subsidence on reclaimed land. When the horizontal displacement is much lower than the vertical displacement, the vertical subsidence dominates the LOS displacements [28]. As XXNA was in flat topography, the displacements acquired from the InSAR method were reasonably interpreted as vertical subsidence in this study. As we can see in Figure 6a, severe subsidence can be observed at XXNA (rectangle A in Figure 6a) and the subsidence rates can reach up to 192 mm/year, which is the position with the largest subsidence rate in the entire study area. In Qiaotou (rectangle B in Figure 6a) and Dengqi (rectangle C in Figure 6a), the maximum subsidence rates can reach 169 mm/year and 156 mm/year, respectively. Except for these three regions mentioned, the rest of Dadeng Island is relatively stable. Figure 6b presents the error map of mean velocity results calculated by the standard deviation value as shown in equation 7. It can be seen that there is high precision in most parts of study area (below ±5 mm/year), and the precision in a very small part of study area (with red color in Figure 6b) is estimated as 5~12 mm/year.
To explore the time series subsidence information on Dadeng Island, five points were selected in Figure 6a. P1, P2, and P3 are located in Region A, while P4 and P5 are distributed in Region B and Region C, respectively. Except for P3 selected on Runway A of XXNA, the other points are undergoing significant subsidence in these regions. Figure 6c presents the time series of the subsidence on these five points. Maximum and minimum cumulative subsidence were observed at P1 and P3, which can reach 213 mm and 55 mm in 486 days, respectively. Furthermore, we can clearly see that P1 to P5 are subsiding linearly. And subsidence trend on P3 is towards stable in the future, while other points will undergo continuous subsidence, which means that the XXNA will probably suffer from severe subsidence problems continuously in the future.  Region C, respectively. Except for P3 selected on Runway A of XXNA, the other points are undergoing significant subsidence in these regions. Figure 6c presents the time series of the subsidence on these five points. Maximum and minimum cumulative subsidence were observed at P1 and P3, which can reach 213 mm and 55 mm in 486 days, respectively. Furthermore, we can clearly see that P1 to P5 are subsiding linearly. And subsidence trend on P3 is towards stable in the future, while other points will undergo continuous subsidence, which means that the XXNA will probably suffer from severe subsidence problems continuously in the future.

Main Subsidence Areas Combined with Land Reclamation and Land Use Planning information
According to the previous studies, ground subsidence is induced by primary consolidation and secondary compression of the clay beneath the reclamation [9,44], and the ground subsidence on Dadeng Island is closely correlated to the completion time of land reclamation [8]. Therefore, we combined historical land reclamation with land use planning to analyze the main subsiding areas on Dadeng Island. The InSAR measurements of regions A, B, and C in Figure 6a are shown in Figures 7-9, respectively. The base map is based on the land use planning of Dadeng Island (see Figure 2). And the dashed lines present boundaries of historical land reclamation (see Figure 3).
in reclaimed land of this region are non-uniform. The maximum subsidence occurred close to the Apron and Runway B in the east of this region with a maximum subsidence rate of 192 mm/year. Combined with historical land reclamation, we found two boundaries of reclaimed land (see the red dashed line and the yellow dashed line in Figure 7) in the region. It indicates that there is a correlation between subsidence and historical land reclamation, i.e., the remarkable subsidence occurred inside the red dashed line and the yellow dashed line, which means that the most significant subsidence was reclaimed in 2014 to 2017. Except for these regions, clear subsidence can also be observed in the west of Runway A, Terminal, Airport Repair Area, and Airport North Road.
As the runways in the airport are the infrastructures that require an extremely stable surface [28], two profiles along the runways were selected to reveal the spatial characteristics of subsidence at XXNA. Profile locations are shown in Figure 7 (i.e., A-A' and B-B'), and the mean subsidence rate along the runways are shown in Figure 8. We can clearly see that there has been obvious significant non-uniform subsidence along profiles A-A' and B-B'. The maximum mean subsidence rate in Runway B can reach 84 mm/year, which is more severe than Runway A (with subsidence of 30 mm/year). It means that the subsidence is uneven along runways and poses a great threat to the safety of aircrafts. To understand the non-uniform subsidence at XXNA, we combined the subsidence with land reclamation history. We marked different land reclamation periods with blue, yellow, and red colors, which indicate the reclamation before 2013, 2013 to 2015, and after 2015. We can see that the subsidence profile curve changes rapidly in the boundary of different reclamation periods. Different reclamation periods mean the land has different physical properties, resulting in different level of subsidence.   historical land reclamation at XXNA. In Qiaotou and Dengqi, the significant subsidence is in the land reclaimed during 2018 to 2019 (see the red dashed line in Figure 9) and during 2016 to 2017 (see magenta dashed line in Figure 10). The significant subsidence occurred in the land which has been reclaimed late in the region of Qiaotou and Dengqi, which confirms that the consolidation and compression in dredger fill are the main cause of subsidence.    Figure 7) in the region. It indicates that there is a correlation between subsidence and historical land reclamation, i.e., the remarkable subsidence occurred inside the red dashed line and the yellow dashed line, which means that the most significant subsidence was reclaimed in 2014 to 2017. Except for these regions, clear subsidence can also be observed in the west of Runway A, Terminal, Airport Repair Area, and Airport North Road.
As the runways in the airport are the infrastructures that require an extremely stable surface [28], two profiles along the runways were selected to reveal the spatial characteristics of subsidence at XXNA. Profile locations are shown in Figure 7 (i.e., A-A' and B-B'), and the mean subsidence rate along the runways are shown in Figure 8. We can clearly see that there has been obvious significant non-uniform subsidence along profiles A-A' and B-B'. The maximum mean subsidence rate in Runway B can reach 84 mm/year, which is more severe than Runway A (with subsidence of 30 mm/year). It means that the subsidence is uneven along runways and poses a great threat to the safety of aircrafts. To understand the non-uniform subsidence at XXNA, we combined the subsidence with land reclamation history. We marked different land reclamation periods with blue, yellow, and red colors, which indicate the reclamation before 2013, 2013 to 2015, and after 2015. We can see that the subsidence profile curve changes rapidly in the boundary of different reclamation periods. Different reclamation periods mean the land has different physical properties, resulting in different level of subsidence. Figures 9 and 10 present detailed spatial subsidence patterns in region B (Qiaotou) and region C (Dengqi), respectively. For Qiaotou area, the most remarkable subsidence reaches 169 mm/year occurred in the northeast of this region, where the land was reclaimed in 2019 and will be used as retail commercial land and residential land. And for Dengqi, the subsidence reaches 156 mm/year in the center region, which was reclaimed in 2017 and will be used as retail commercial land, commercial and financial land. It should be noted that there will be a metro built in Qiaotou and Dengqi. Slight subsidence was observed near Dadeng Station in Qiaotou (see Figure 9), while the Shuanghu Station in Dengqi (see Figure 10) has significant subsidence. In addition, the observed subsidence in Qiaotou will affect the Xi East Road and Yang Tang Road (see solid white lines in Figure 9) and in Dengqi, the subsidence will affect Yingbing Road and Airport Hightway (see Figure 10). The reclaimed boundaries in Qiaotou and Dengqi also show a correlation between subsidence and historical land reclamation at XXNA. In Qiaotou and Dengqi, the significant subsidence is in the land reclaimed during 2018 to 2019 (see the red dashed line in Figure 9) and during 2016 to 2017 (see magenta dashed line in Figure 10). The significant subsidence occurred in the land which has been reclaimed late in the region of Qiaotou and Dengqi, which confirms that the consolidation and compression in dredger fill are the main cause of subsidence.

Discussion
To quantitatively analyze the potential risk caused by subsidence in the land use planning, we combined subsidence monitoring results with twenty-two land use types according to future land use planning (XXNA was divided into 5 different land use types). According to four subsidence rate classes and the Expectation and Entropy methods, we calculated the area and the values of Ex and En in different land use types. The statistical results are shown in Table 2. We divided different land use types of Dadeng Island into four risk classes by subsidence rate to show the spatial distribution of potential risk caused by subsidence on Dadeng Island (Figure 11a). High risk areas, accounting for about 3.87% of the total land area (Figure 11b), are mainly distributed in Land for Airport, Road, Retail Commercial Land, Apron, Green Space and Square Land, Class I and II Residential Land, and Commercial and Financial Land (see Table 2). Moderate risk areas account for 2.42% of all and mainly in Land for Airport, Road, Green Space and Square Land, and Class I and II Residential Land (see Table 2). Low risk areas and Very Low risk areas represent about 8.32% and 85.39%, which means that most areas of Dadeng Island are stable and suitable for future land use.   Figure 6a).

Discussion
To quantitatively analyze the potential risk caused by subsidence in the land use planning, we combined subsidence monitoring results with twenty-two land use types according to future land use planning (XXNA was divided into 5 different land use types). According to four subsidence rate classes and the Expectation and Entropy methods, we calculated the area and the values of Ex and En in different land use types. The statistical results are shown in Table 2. We divided different land use types of Dadeng Island into four risk classes by subsidence rate to show the spatial distribution of potential risk caused by subsidence on Dadeng Island (Figure 11a). High risk areas, accounting for about 3.87% of the total land area (Figure 11b), are mainly distributed in Land for Airport, Road, Retail Commercial Land, Apron, Green Space and Square Land, Class I and II Residential Land, and Commercial and Financial Land (see Table 2). Moderate risk areas account for 2.42% of all and mainly in Land for Airport, Road, Green Space and Square Land, and Class I and II Residential Land (see Table 2). Low risk areas and Very Low risk areas represent about 8.32% and 85.39%, which means that most areas of Dadeng Island are stable and suitable for future land use.
As the areas of different land use types have significantly different subsidence rate, Figure 12 presents the percentage of potential risk in different land use types. From the results, we can clearly see that the potential risk in Special Land is highest, which has 13% high risk and 26% moderate risk. Then there is Sports Land, which has 13% high risk and 17% moderate risk. In addition, the high risk areas in Commercial and Financial Land, Apron, and Retail Commercial Land are more than 5%. For these land use types, subsidence would have significant influence to the future land use and the potential subsidence risk should be considered. In contrast, the land use types of Other Public Management and Public Service Facilities and Class III Residential Land have very low risk, which are suitable for future land use. For XXNA, the highest potential risk is in Apron and Airport Infrastructure, which have 8% high risk and 4% moderate risk, and 5% high risk and 7% moderate risk, respectively. The moderate potential risk is in Land for Airport and Runway, which have 3% high risk and 2% moderate risk, and 2% high risk and 3% moderate risk, respectively. The lowest potential risk is in Terminal, which only has 6% low risk and 94% very low risk. Therefore, attention has to be paid to reduce the potential harm of the subsidence in future land use of Dadeng Island, especially for XXNA, Special Land and Sports Land.  As the areas of different land use types have significantly different subsidence rate, Figure 12 presents the percentage of potential risk in different land use types. From the results, we can clearly see that the potential risk in Special Land is highest, which has 13% high risk and 26% moderate risk. Then there is Sports Land, which has 13% high risk and 17% moderate risk. In addition, the high risk areas in Commercial and Financial Land, Apron, and Retail Commercial Land are more than 5%. For these land use types, subsidence would have significant influence to the future land use and the potential subsidence risk should be considered. In contrast, the land use types of Other Public Infrastructure, which have 8% high risk and 4% moderate risk, and 5% high risk and 7% moderate risk, respectively. The moderate potential risk is in Land for Airport and Runway, which have 3% high risk and 2% moderate risk, and 2% high risk and 3% moderate risk, respectively. The lowest potential risk is in Terminal, which only has 6% low risk and 94% very low risk. Therefore, attention has to be paid to reduce the potential harm of the subsidence in future land use of Dadeng Island, especially for XXNA, Special Land and Sports Land.  Figure 13 presents Expectation (Ex) and Entropy (En) values in different land use types. As mentioned in Section 3.2, Ex and En values can be used to analyze the overall level and the degree of non-uniformity land subsidence in terms of different land use. The results show significant difference of potential risk. Some land use types, such as Special Land and Sports Land, have high Ex values (more than -30) and high En values (more than 20). It means that the subsidence magnitude and the non-uniformity degree of subsidence are large in these areas. While for some other land use types, such as Regional Transportation Facilities Land, have low Ex values (Ex value below 5) and moderate En values (more than 9), i.e., the subsidence magnitude is low, but the non-uniform degree is large. The non-uniform subsidence is the most serious problem in these areas. Other Public Management and Public Service Facilities and Class III Residential Land have low values in Ex and En (Ex value below 5 and En value below 3). These areas are quite stable and have no potential subsidence risk. Therefore, the overall level and the non-uniformity degree of land subsidence in different land use types are determined by Ex and En methods, which provide a qualitative reference for prevention and control of subsidence geohazard.  As mentioned in Section 3.2, Ex and En values can be used to analyze the overall level and the degree of non-uniformity land subsidence in terms of different land use. The results show significant difference of potential risk. Some land use types, such as Special Land and Sports Land, have high Ex values (more than -30) and high En values (more than 20). It means that the subsidence magnitude and the non-uniformity degree of subsidence are large in these areas. While for some other land use types, such as Regional Transportation Facilities Land, have low Ex values (Ex value below 5) and moderate En values (more than 9), i.e., the subsidence magnitude is low, but the non-uniform degree is large. The non-uniform subsidence is the most serious problem in these areas. Other Public Management and Public Service Facilities and Class III Residential Land have low values in Ex and En (Ex value below 5 and En value below 3). These areas are quite stable and have no potential subsidence risk. Therefore, the overall level and the non-uniformity degree of land subsidence in different land use types are determined by Ex and En methods, which provide a qualitative reference for prevention and control of subsidence geohazard.

Conclusions
In this paper, the time series InSAR technique is used to monitor reclaimed land subsidence at XXNA for obtaining the spatiotemporal evolution information of subsidence. The InSAR measurement results show that some of the reclaimed areas at XXNA are still in high subsidence rates with a maximum of -192 mm/year from January 2018 to April 2019. Three main subsidence areas are identified, i.e., runway region, Qiaotou, and Dengqi area. Specially, the high linearly-rate subsidence in the main subsidence area in this paper is consistent with previous studies, which means those areas will probably undergo continuous subsidence in the future. These findings prove that the SAR interferometry is a powerful tool for evaluating the land subsidence geohazards in the reclaimed land.
By analyzing the spatial subsidence characterizations of the main subsiding areas combined with historical land reclamation and future land use planning, we find that the significant non-uniform subsidence poses a potential threat to the safety of the future land use. In addition, clear correlations between historical reclamation activities and land subsidence have been discovered, which confirms that the consolidation and compression in dredger fill is the main cause of subsidence. Regional land subsidence and its differences will directly affect future land use on reclaimed land. Therefore, we combined subsidence monitoring results with twenty-two land use types in land use planning and used Ex and En methods to characterize potential risk caused by subsidence in study area. Results showed that 3.4% and 2.4% of land were suffering from high and moderate subsidence risk, respectively. Different land use types showed significant difference in potential risk. In future work, as the subsidence is likely to continue, it is suggested to perform continuous monitoring of the XXNA.

Conclusions
In this paper, the time series InSAR technique is used to monitor reclaimed land subsidence at XXNA for obtaining the spatiotemporal evolution information of subsidence. The InSAR measurement results show that some of the reclaimed areas at XXNA are still in high subsidence rates with a maximum of -192 mm/year from January 2018 to April 2019. Three main subsidence areas are identified, i.e., runway region, Qiaotou, and Dengqi area. Specially, the high linearly-rate subsidence in the main subsidence area in this paper is consistent with previous studies, which means those areas will probably undergo continuous subsidence in the future. These findings prove that the SAR interferometry is a powerful tool for evaluating the land subsidence geohazards in the reclaimed land.
By analyzing the spatial subsidence characterizations of the main subsiding areas combined with historical land reclamation and future land use planning, we find that the significant non-uniform subsidence poses a potential threat to the safety of the future land use. In addition, clear correlations between historical reclamation activities and land subsidence have been discovered, which confirms that the consolidation and compression in dredger fill is the main cause of subsidence. Regional land subsidence and its differences will directly affect future land use on reclaimed land. Therefore, we combined subsidence monitoring results with twenty-two land use types in land use planning and used Ex and En methods to characterize potential risk caused by subsidence in study area. Results showed that 3.4% and 2.4% of land were suffering from high and moderate subsidence risk, respectively. Different land use types showed significant difference in potential risk. In future work, as the subsidence is likely to continue, it is suggested to perform continuous monitoring of the XXNA. Based on the combination of key subsidence areas and the land use types, Road, Apron, Special Land, and Sports Land should receive special attention with detailed monitoring. The risk prevention plan should be considered in advice for the areas with potential subsidence geohazards. This work gives a guideline for the future ground subsidence monitoring at XXNA and provides a reference for potential hazard evaluation based on the land use design and high subsidence rate area.