Diagnosing Subsidence Geohazard at Beijing Capital International Airport, from High-Resolution SAR Interferometry

Beijing Capital International Airport (BCIA) has suffered from uneven land subsidence since 1935, which affects the smoothness of airport runways and seriously threatens the safety of aircrafts. In this paper, a spaceborne interferometric synthetic aperture radar (InSAR) with high-resolution Cosmo-SkyMed SAR data was utilized at BCIA for the first time to diagnose the subsidence hazard. The results show that subsidence is progressing at BCIA at a maximum rate of 50 mm/year, which is mainly distributed in the northwest side of the airport. It was found that the Shunyi-Liangxiang fault directly traverses Runway2 and Runway3 and causes uneven subsidence, controlling the spatial subsidence pattern to some degree. Four driving factors of subsidence were investigated, namely: the over-exploitation of groundwater, active faults, compressible soil thickness, and aquifer types. For the future sustainable development of BCIA, the influence of Beijing new airport and Beijing Daxing International Airport (BDIA), was analyzed and predicted. It is necessary to take relevant measures to control the uneven subsidence during the initial operation of BDIA and conduct long-term monitoring to ensure the regular safe operation of BCIA. This case demonstrates a remote sensing method of diagnosing the subsidence hazard with high accuracy and non-contact, providing a reliable alternative for the geohazard diagnosis of key infrastructures in the future.


Introduction
Beijing, the political center of China with a population of more than 21 million, has suffered from land subsidence since 1935 [1]. Severe subsidence still occurs in eastern Beijing, with a maximum subsidence rate exceeding 10 cm/year [2][3][4]. Beijing Capital International Airport (BCIA) is located in the Beijing Shunyi District, northeast of downtown, which was involved in the core subsidence area in Beijing. The runways in the airport are the infrastructures that most require stable ground [5].

Study Area, Data Source
The Beijing Plain is located in the northern part of the North China Plain, which is formed by sediments carried by Chaobai River and Wenyu River, surrounded by the Yanshan Mountains in the north and the Taihang Mountains in the west. It is a typical piedmont alluvial plain and most of the area is covered by fine sediments (Figure 1a). BCIA, located 32 km northeast of Beijing's city center, is the 4F-level international airport serving Beijing. The airport is located in a flat zone between the Wenyu River and the Chaobai River alluvial fan at the end of the Yongding River alluvial fan, with an average elevation of 30.4 m (Figure 1a). It has been the world's second busiest airport in terms of passenger traffic since 2010. In 2018, the airport handled of 614-thousand air traffic movements, 100 million passengers, and over 2 million tons of cargo in total [24], becoming the second airport in the world with an annual passenger throughput of over 100 million.
BCIA was opened in 1958 and three round expansion were performed in 1980, 1999 and 2007. Now there are in total three terminals and three runways within the airport (Figure 1b). Terminal 1 (T1) and Terminal 2 (T2), in the west part of the airport linked by a public walkway, were opened on 1980 and 1999 with a floor area 60000 m 2 and 336000 m 2 [25], respectively. Terminal 3 (T3), opened in 2008, is the largest airport terminal building complex built in a single phase, with 986,000 m 2 in total floor area, consisting of a main passenger terminal (Terminal 3C, T3C), two satellite concourses (Terminal 3D and Terminal 3E, T3D and T3E), and five floors above ground and two underground [25] ( Figure 1b). There are, in total, three runways in the BCIA. All of them have the same geographic bearing and magnetic bearing of 173/353 • and 179/359 • , respectively [24]. The dimensions of Runway2 and Runway3 are 3800 × 60 m and 3200 × 50 m for Runway1.  (Terminal 3D and Terminal 3E, T3D and T3E), and five floors above ground and two underground [25] (Figure 1b). There are, in total, three runways in the BCIA. All of them have the same geographic bearing and magnetic bearing of 173/353° and 179/359°, respectively [24]. The dimensions of Runway2 and Runway3 are 3800 × 60 m and 3200 × 50 m for Runway1.
CSK is an Earth-observation satellite space-based radar system funded by the Italian and Ministry of Defence (MoD), controlled and maintained by the Italian Space Agency (ASI), intending for Dual-use (both military and civilian use) to establish a global service supplying provision of SAR data [26]. The space segment of the system includes a constellation of four mid-sized satellites and each of them was equipped with a multi-mode high-resolution SAR operating at X-band with global coverage [27]. Each satellite repeats the same ground track every 16 days, and all of the satellites follow the same ground track of between 1 and 15 days [26], providing observation of an area of interest with 1-16 repeated cycle in all-weather conditions. The CSK satellites have three basic types of imaging modes, namely Spotlight mode, Stripmap mode, ScanSAR mode.
In this study, twenty-one CSK images were acquired during November 2017 to March 2019 with descending orbit in Stripmap mode. The images are single polarization (HH) data in SLC (single look complex) format, with a spatial resolution of 3 m and an incident angle of 20.1°. The wavelength length is 3.1 cm, which could measure small displacement for InSAR process. The dates, temporal baseline and polarization mode of the acquisitions were listed in the Table 1, while the spatial perpendicular baselines of 21 CSK acquisitions are shown in Figure 2. CSK is an Earth-observation satellite space-based radar system funded by the Italian and Ministry of Defence (MoD), controlled and maintained by the Italian Space Agency (ASI), intending for Dual-use (both military and civilian use) to establish a global service supplying provision of SAR data [26]. The space segment of the system includes a constellation of four mid-sized satellites and each of them was equipped with a multi-mode high-resolution SAR operating at X-band with global coverage [27]. Each satellite repeats the same ground track every 16 days, and all of the satellites follow the same ground track of between 1 and 15 days [26], providing observation of an area of interest with 1-16 repeated cycle in all-weather conditions. The CSK satellites have three basic types of imaging modes, namely Spotlight mode, Stripmap mode, ScanSAR mode.
In this study, twenty-one CSK images were acquired during November 2017 to March 2019 with descending orbit in Stripmap mode. The images are single polarization (HH) data in SLC (single look complex) format, with a spatial resolution of 3 m and an incident angle of 20.1 • . The wavelength length is 3.1 cm, which could measure small displacement for InSAR process. The dates, temporal baseline and polarization mode of the acquisitions were listed in the Table 1, while the spatial perpendicular baselines of 21 CSK acquisitions are shown in Figure 2.

Methodology
Small Baseline Subsets InSAR (SBAS-InSAR) technology utilizes a number of SAR images covering the same area to extract the ground target points with high coherence. Based on the differential phase of these points and the subsequent modelling, the precise time series displacements information can be retrieved. The whole process can be divided into 4 parts, as shown in the flowchart in Figure 3.

Methodology
Small Baseline Subsets InSAR (SBAS-InSAR) technology utilizes a number of SAR images covering the same area to extract the ground target points with high coherence. Based on the differential phase of these points and the subsequent modelling, the precise time series displacements information can be retrieved. The whole process can be divided into 4 parts, as shown in the flowchart in Figure 3. Assuming that the N+1 SAR images covering the same area were obtained at the ordered time .., n t t t + , in the data pre-processing, the raw SLC SAR data can be cropped and coregistrated into target SLC datasets with the precise orbit. In the interferometric process, M interferometric pairs Assuming that the N+1 SAR images covering the same area were obtained at the ordered time (t 0 , t 1 , . . . , t n+1 ), in the data pre-processing, the raw SLC SAR data can be cropped and coregistrated into target SLC datasets with the precise orbit. In the interferometric process, M interferometric pairs can be selected from the total N (N+1)/2 combination based on the given temporal and spatial threshold values to form a baseline network. Then the interferometric processing includes interferogram generation, flat and topographic phase removal and adaptive filter, in which the external DEM will be involved. Based on the 3D phase unwrapping algorithm, the phase unwrapping for the M interferogram is performed in order to obtain the unwrapped phase change in the line of sight (LOS) direction.
The unwrapped phase on the arbitrary pixels (x, y) in the j th interferogram generated at times t A and t B can be modeled as [28], atm (x, r) denotes the displacement in LOS direction, topographic error caused by the inaccuracy of external DEM and the atmospheric artifacts, respectively. φ j de f (x, r) could be expressed as the product of the mean phase velocity between time adjacent acquisitions multiplying time span, i.e., where λ, B ⊥ , r, θ is the radar wavelength, perpendicular baseline, range distance and incidence angle, respectively. Therefore, by combining Equations (1) to (3), the first time series estimation can be performed by, N unknown deformation rates and the topographic error Z(x, r) can be solved in a least squares sense (if all the acquisitions are well connected, we should have an overdetermined system, i.e., M > N+1, or the singular value decomposition (SVD) is performed). Then the atmospheric components can be estimated and eliminated using filters through the residue phase after first estimation. In the final estimation using Equation (4), the refined time series displacements, topographic correction and average displacement rate results would be obtained using a least-squares method in the time domain. Finally, the displacement results were geocoded into the WGS84 UTM coordinate system, projected to subsidence and exported for post-processing analysis. The details of the SBAS-InSAR method could be found from [28][29][30][31].

Results
Based on the high-resolution CSK images and the SBAS-InSAR technique, the subsidence of BCIA during November 2017 to March 2019 was acquired. According to the previous studies (e.g., [2]), the horizontal displacement was less than 1 mm/year in Beijing, which is much lower than the vertical displacement. It means that the vertical subsidence dominates the LOS displacements. Therefore, as BCIA was in flat topography, the displacements acquired from the InSAR technique were reasonably interpreted as vertical subsidence in this study. A mean velocity subsidence map was acquired (Figure 4), in which the negative value indicates the targets moving far away from the satellites (i.e., subsidence). Uneven spatial subsidence distribution can be observed. Severe subsidence was mainly distributed in the north and west of BCIA, while the middle-east and south-east sections of BCIA were relative stable. The severe subsidence affected the group of buildings in the north of T1 and T2, and the north part of Runway2 with a maximum subsidence rate of 50 mm/year. In terms of the terminals, the subsidence on T1 and T2 was about 35 mm/year, while T3 was relative stable. It should be noted that although the severe subsidence occurred in the north of Runway2, the middle and south part of Runway2 are stable, indicating that extensive uneven subsidence existed on this runway. Runway3 is stable, which is seen from its high-coherence points-which look to be much fewer in quantity than those on and Runway1 and Runway2. The detailed subsidence along runways and terminals will be presented in the following paragraphs.  In order to reveal the uneven subsidence affecting BCIA, the subsidence on T2 and Runway2 was enlarged; shown in Figures 5a and 5b, respectively. When interpreting the color high-coherence points in line with subsidence information, it is clear that the uneven subsidence occurred on T2. The north and middle sections of T2 were suffering from subsidence at a rate of more than 30 mm/year, while the rate for the south section was 17 mm/year. A bigger subsidence occurred on Runway2. The subsidence in the north of Runway2 reached up to 44 mm/year, while it was 11 mm/year in the middle and south of Runway2. It should be noted that the subsidence difference was more than 30 mm/year in the are 2000 m north of the middle section of Runway2, which would pose a potential risk to the safety of the aircraft during take off and landing. In order to reveal the uneven subsidence affecting BCIA, the subsidence on T2 and Runway2 was enlarged; shown in Figure 5a,b, respectively. When interpreting the color high-coherence points in line with subsidence information, it is clear that the uneven subsidence occurred on T2. The north and middle sections of T2 were suffering from subsidence at a rate of more than 30 mm/year, while the rate for the south section was 17 mm/year. A bigger subsidence occurred on Runway2. The subsidence in the north of Runway2 reached up to 44 mm/year, while it was 11 mm/year in the middle and south of Runway2. It should be noted that the subsidence difference was more than 30 mm/year in the are 2000 m north of the middle section of Runway2, which would pose a potential risk to the safety of the aircraft during take off and landing. The time series results on some selected typical points were shown in Figure 6. Points #1, #2 and #5 are in the main subsidence area of BCIA with an accumulated subsidence of more than 60 mm from November 2017 to April 2019. The points #3 and #4 are in the other side of T2 and Runway2 with a relative minor accumulated subsidence of less than 30 mm and 20 mm, respectively, indicating the uneven subsidence on T2 and Runway2. The points #6 and #7 were on the T3 and Runway3, where the subsidence is relatively small, with an accumulated subsidence of less than 10 mm. From the time series results it can be seen that the subsidence was slightly accelerated from May 2018 to August 2018, possibly as a result of the rainy season in Beijing during this time period.
It can be observed from Figure 4 that a linear boundary can be found in the subsidence spatial pattern. After the subsidence results were interpolated, as shown Figure 7a, it is obvious that at two sides of this linear boundary, a clear subsidence difference can be recognized. This linear boundary spatially corresponds to the Shunyi-Liangxiang fault trace (the projections of fault onto the ground surface), indicating a potential link between them. The subsidence in the northwest of the fault trace is much more severe than it in the southeast of the fault trace. It should be noted this fault trace and linear boundary traverses the whole airport, including the Runway2 and Runway3, already causing uneven subsidence for the runways. The time series results on some selected typical points were shown in Figure 6. Points #1, #2 and #5 are in the main subsidence area of BCIA with an accumulated subsidence of more than 60 mm from November 2017 to April 2019. The points #3 and #4 are in the other side of T2 and Runway2 with a relative minor accumulated subsidence of less than 30 mm and 20 mm, respectively, indicating the uneven subsidence on T2 and Runway2. The points #6 and #7 were on the T3 and Runway3, where the subsidence is relatively small, with an accumulated subsidence of less than 10 mm. From the time series results it can be seen that the subsidence was slightly accelerated from May 2018 to August 2018, possibly as a result of the rainy season in Beijing during this time period.
It can be observed from Figure 4 that a linear boundary can be found in the subsidence spatial pattern. After the subsidence results were interpolated, as shown Figure 7a, it is obvious that at two sides of this linear boundary, a clear subsidence difference can be recognized. This linear boundary spatially corresponds to the Shunyi-Liangxiang fault trace (the projections of fault onto the ground surface), indicating a potential link between them. The subsidence in the northwest of the fault trace is much more severe than it in the southeast of the fault trace. It should be noted this fault trace and linear boundary traverses the whole airport, including the Runway2 and Runway3, already causing uneven subsidence for the runways. Subsidence profiles along the three runways were presented in Figure 7b. The profiles AA', BB' and CC' represent the subsidence along Runway1, Runway2 and Runway3, respectively. It can be observed that Runway1 (profile AA') was near the subsidence core area with a maximum subsidence rate of up to 40 mm/year. The subsidence grows more severe from north to south, i.e., the subsidence increased from 20 mm/year to 38 mm/year. The subsidence along Runway1 was not affected by the fault. The Shunyi-liangxiang fault traverses the whole Runway2 and cut it into two part. The location of Shunyi-liangxiang fault was denoted by double vertical line in Figure 7b, where a dramatic change of subsidence was observed. It is indicated that the Shunyi-liangxiang fault has significant influence on the subsidence and caused −40 mm/year uneven subsidence to Runway2. The Shunli-Liangxiang Fault is a normal fault, which may be useful for explaining how the subsidence trend corresponds to expected compressional or extensional deformation on the fault. There is no clear uneven subsidence along Runway3, although the fault also traversed a corner of Runway3. Runway3 was relatively stable with a maximum subsidence rate of 10 mm/year, which occurred in the south end. These results show a good agreement with the results derived in 2003-2017 (i.e., [2,23,24]). Compared with these studies, the subsidence rate is slower, indicating that valid relevant measures have been taken to slow the subsidence in 2017-2019. Subsidence profiles along the three runways were presented in Figure 7b. The profiles AA', BB' and CC' represent the subsidence along Runway1, Runway2 and Runway3, respectively. It can be observed that Runway1 (profile AA') was near the subsidence core area with a maximum subsidence rate of up to 40 mm/year. The subsidence grows more severe from north to south, i.e., the subsidence increased from 20 mm/year to 38 mm/year. The subsidence along Runway1 was not affected by the fault. The Shunyi-liangxiang fault traverses the whole Runway2 and cut it into two part. The location of Shunyi-liangxiang fault was denoted by double vertical line in Figure 7b, where a dramatic change of subsidence was observed. It is indicated that the Shunyi-liangxiang fault has significant influence on the subsidence and caused −40 mm/year uneven subsidence to Runway2. The Shunli-Liangxiang Fault is a normal fault, which may be useful for explaining how the subsidence trend corresponds to expected compressional or extensional deformation on the fault. There is no clear uneven subsidence along Runway3, although the fault also traversed a corner of Runway3. Runway3 was relatively stable with a maximum subsidence rate of 10 mm/year, which occurred in the south end. These results show a good agreement with the results derived in 2003-2017 (i.e., [2,23,24]). Compared with these studies, the subsidence rate is slower, indicating that valid relevant measures have been taken to slow the subsidence in 2017-2019.

Driving Factors of Subsidence
The stability of ground is very important for the regular operation of airport. Uneven subsidence would bring potential risk to the runways and terminals, especially for the aircraft during the take-

Driving Factors of Subsidence
The stability of ground is very important for the regular operation of airport. Uneven subsidence would bring potential risk to the runways and terminals, especially for the aircraft during the take-off and landing. Uneven subsidence has occurred in dozens of airports built on land reclaimed from the sea, for example, Lap Kok Airport (Hong Kong, China) [32], Kansai International Airport (Osaka, Japan) [33,34], San Francisco International Airport (San Francisco, USA) [35]. Based on the previous study, the over extraction of groundwater is the main cause for the large-scale land subsidence in Beijing [8,36], resulting in BCIA being a very typical inland airport suffering from subsidence at a maximum rate of 50 mm/year.
Beijing has been suffering from land subsidence since 1935 and is expected to continue to. The severe subsidence areas in Beijing are mainly distributed over eastern Beijing, i.e., the east of Chaoyang District and the north of Tongzhou District with a maximum rate up to 140 mm/year till 2017 [2]. BCIA is in an enclave of Chaoyang District and the surroundings of that enclave in suburban Shunyi District, which is 10 km away to the severe subsidence area. Therefore, the wide-scale subsidence in Beijing, investigated in [4,8,37], will definitely have an effect on the subsidence in BCIA. Some driving factors causing the subsidence were analyzed and listed as following and all of them should be taken into account to take measures to slow down the subsidence.

Over-Exploitation of Groundwater
The aquifer system is composed by aquifers and the aquitards. When groundwater is pumped from the underground aquifer system, it undergoes deformation to different degrees. Most of the permanent subsidence occurs due to the irreversible compression or consolidation of aquitards [38]. Thus, the deformation of an aquifer-system can be inelastic or elastic. The over-exploitation of groundwater is supposed as the main cause of land subsidence [36]. In terms of time series comparison between subsidence and the underground water level, a similar elastic behavior trend can be found [38]. Therefore, controlling the exploitation of groundwater will play a dominant role in controlling subsidence in Beijing.

Active Faults
Beijing is located in the faults junction region between the faults system in Yanshan, North China Plain and Taihang Mountain [39]. The fault systems mainly consist of the NE faults system (e.g., Nankou Piedmont fault, XiaotangMountain-Dongbeiwang fault, Huangzhuang-Gaoliying fault) and the NW fault systems (e.g., NankouSunhe and Yongding River faults). As revealed by the above results and previous studies [2], it is clear that the subsidence pattern is spatially controlled by geological faults to some degree. In particular, the Shunyi-liangxiang fault traverses the whole BCIA, causing the non-negligible uneven subsidence over Runway2 and T2.

Compressible Soil Thickness
In general, when piezometric levels decrease, the thicker the accumulated compressible soil is, the greater the land subsidence rates are [40]. It was revealed that in same underground water exploitation conditions, the clayey units (aquitards made of fine soils) are two or three times more deformable than those of sand and gravel (coarse soils) [41]. Figure 8 presents the statistics on the distribution of subsidence in different thicknesses of soft soil, i.e., the percentage of each compressible thickness in different subsidence rate. The subsidence in Figures 8 and 9 was derived by InSAR from 2003 to 2010, and the statistics were based on [8]. It is clear that when the thickness of soft soil is over 80 m, the subsidence is less than −40 mm/year. The BCIA is in the soft soil with a thickness of 60-70 m, which means that the potential severe subsidence up to 100 should be paid attention.

Aquifer Types
In the Beijing plain, the structure of the aquifer system is composed of three layers, namely coarse grain sediments with single sandy gravel, median grain sediments with multiple layers of pebbles and gravels, and fine grain sediments with multiple sand layers and gravel layers [36]. The percentage of each grain sediment in the different subsidence rates is shown in Figure 9. It is obvious that the subsidence over −40 mm/year almost exclusively occurred in the fine grain sediment area. Unfortunately, the whole northeastern Beijing are in fine grain sediment, which provide potential risk for the subsidence developing.

Future Sustainable Development
With rapid economic growth and urbanization, the huge population in Beijing has placed further stress on the existing transport facilities, including air travel. Over the next few years, a large number of passengers are expected to travel in and out of the country through the airport. The BCIA is currently running at maximum capacity and cannot be expanded further [42]. The growth of the city has necessitated the construction of a new airport.
Beixing Daxing International Airport (BDIA) was approved for construction in January 2013, situated in Daxing that a suburb district 46 kilometers south from the city center of Beijing. The new airport is incurred an approximate expenditure of CNY 120 billion (USD 17.46 billion) [42]. BDIA is the first terminal in the world that adopted a two-story design for departures to address overflowing aviation traffic congestion, as an international hub connecting to the world and a new engine for China's economy. Beijing Nanyuan Airport (BNA) is the first airport in Chinese history located around 13 km from the city center in Fengtai District. It was put into civil use in 2007 after the first reconstruction of the airport. BNA is a relatively small airport serving Beijing and will be closed following the commencement of the full operation of BDIA.
It is certain that BDIA will have influence on the throughput of BCIA. Figure 10 presents the current and predicted total number of flights from three different airports in Beijing, i.e., BCIA, BDIA and BNA. It can be seen that the civil flights will not depart from BNA again once the official operation of BDIA has begun. After that, the number of flights in BCIA will decrease gradually as those leaving and entering BDIA increase dramatically. By 2021, the number of flights in BCIA and BDIA will reach up to the same level with a difference of less than three hundred flights. Compared with the period before the operation of BDIA, the flights in BCIA in 2021 should decrease by more than four hundred (29.63%). The current and predicted total throughputs of the different airports in Beijing are presented in Figure 11. In 2020, with the official operation of BDIA, the throughput in BCIA will dramatically decrease by 17 million, while there will be a 25 million increase in the throughput at BDIA. According to the prediction, the total throughput of BDIA will rapidly increase to 72 million in the 6 years from 2019 to 2025. The throughput of BCIA will decrease to 66 million in 2021, and a rebound is expected from 2021 to 2025, finally reaching up to 80 million in 2025.     Based on the above analysis, the total number of flights and throughput in BCIA will be significantly influenced by BDIA. BCIA has an annual passenger handling capacity of 75 million [42]. In 2018, BCIA handled over 100 million passengers, exceeding its capacity by 33% and reaching full saturation. The operation of BDIA will reduce congestion, the stress and load on the existing airport facilities in BCIA, e.g., the terminals and runways. It should be noted that the throughput at BCIA will again trend upwards after 2021. Therefore, the related restoration work on the subsidence of terminals and runways in BCIA should begin work in 2020-2021.

Conclusions
In this paper, the current subsidence at BCIA from November 2017 to April 2019 was revealed using satellite SAR interferometry based on 21 CSK images. The results show that the subsidence on BCIA continued in 2017-2019, with a maximum subsidence rate of more than 50 mm/year. Compared with previous studies, the rate has slowed down to some degree. The high-resolution CSK data enables the detailed analysis of the subsidence on the terminals and runways, and its relationship with geological faults. It was found that the uneven subsidence extensively affected the Terminal 2 and Runway2 at BCIA. More importantly, in the spatial subsidence map, it can be seen that the subsidence pattern exhibits two parts that are separated by a near-linear boundary, which corresponds to the Shunyi-Liangxiang fault. A clear difference in subsidence can be observed at the two sides of this fault, indicating that the fault may have spatially controlled the subsidence at BCIA to some degree. It should be noted that the profile along Runway2 indicates that it was significantly affected by the uneven subsidence, with a different subsidence rate of 30 mm/year within 2 km.
To analyze the future sustainable development of BCIA, four driving factors affecting the subsidence were investigated, i.e., the over-exploitation of groundwater, active faults, compressible soil thickness and aquifer types. The influence of BDIA was analyzed by comparing the current and predicted total number of flights and the throughput of three airports in Beijing. It is necessary to take relevant measures to control the uneven subsidence during the operation of BDIA, and to conduct long-term monitoring at BCIA to ensure the safety of aircraft and safe regular operation. This case demonstrates a remote sensing method of diagnosing the subsidence hazard, with a high accuracy and non-contact for Geohazard and risk reduction in the future.