Quantitative Evaluations of Pumping-Induced Land Subsidence and Mitigation Strategies by Integrated Remote Sensing and Site-Specific Hydrogeological Observations

1. Introduction

Land subsidence, the gradual settlement of the ground surface on a large scale, is an environmental issue caused by various natural or anthropogenic factors [1,2,3]. Tectonic activities, soil consolidation, or other physicochemical processes might trigger land subsidence naturally. On the other hand, human activities, including groundwater or oil extraction, mining, or construction loading, are also the causes of ground surface settlement. Among these causes, long-term groundwater extraction stands out as the primary trigger of land subsidence, especially in populated and developed areas where the geological foundation is constituted by the recent alluvial, marine, or lacustrine deposits [2,4,5]. The consequences of land subsidence are serious since it significantly affects the buildings or transport infrastructures and might lead to inundation in low-lying areas, resulting in heavy financial burdens for local authorities. Cities or regions like the San Joaquin Valley (San Francisco) [6,7], Mexico City (Mexico) [8], Beijing (China) [9], Hanoi (Vietnam) [10], and Pingtung County (Taiwan) [11], among others, are examples of the land subsidence due to excessive groundwater extraction.

The Choushui River Fluvial Plain (CRFP), one of the main agricultural areas in Taiwan, is no exception. This region is well-known for severe land subsidence due to massive groundwater exploitation [12,13,14,15], which has been an urgent issue since it might potentially pose a serious threat to the Taiwan High-speed Rail (THSR) passing through this region [16]. The groundwater exploitation leads to sediment compression [17], which might trigger an angular deflection if it occurs below the foundation of the railway pillars [18,19,20], severely affecting the THSR operation safety. Because the THSR is a continuous railway system, the breakdown at any section might stagnate the entire operation; hence, THSR safety standards must be a top priority [16]. Therefore, establishing a monitoring network is necessary, especially near the THSR, to promptly obtain the subsidence progress, based on which the local authorities can introduce appropriate policies to control the influencing factors and mitigate the sinking rates.

The subsidence monitoring network previously utilized traditional geodetic methods, such as the Continuous Global Positioning System (CGPS) and precise leveling survey. However, the sparse density of measuring points and the time between sampling constrain the effectiveness of these techniques [21]. For example, the measurements provided by these methods are often interpolated to generate a deformation map. The interpolation process might introduce errors proportional to the distance among stations, especially when missing data exists [22]. The Differential Interferometric Synthetic Aperture Radar (DInSAR) techniques provide alternatives to overcome the intrinsic limitations. This DInSAR technique extracts the surface deformation in the phase difference of at least two SAR images acquired on different days over an identical area [23] and provides wide coverage measurements. However, the DInSAR performance is significantly impaired by temporal decorrelation and atmospheric delay [24]. It is recognized that the paddy fields, whose surface properties vary seasonally, primarily cover the land surface of the CRFP and may cause spatiotemporal decorrelation effects, as previously mentioned [13,25]. For that reason, the time-series InSAR techniques, including Permanent Scatter Interferometry (PSI) [26,27] and Small Baseline Subset (SBAS) [28], which exploit a great number of SAR images acquired in the same area over time and extract information from pixels stable in scattering properties, are developed to reduce the effects of such degradation problems [23,29].

For decades, several studies have estimated the subsidence rate in the CRFP and provided the subsidence profiles along the THSR by applying either PSI or SBAS method. Hung, et al. [14] exploited 20 ENVISAT SAR images, then integrated the PSI results and leveling measurements to estimate the vertical displacements during 2006–2008, which showed that three townships of Tuku, whose annual subsidence rates were up to 70 mm/year in maximum, were under the railway of THSR. Lu, et al. [30] combined PSI results extracted from 34 ENVISAT SAR images using geostatistical methods and validated them with GPS data. Then, the earlier-derived deformation patterns were compared with three aquifers’ annual groundwater level fluctuations to analyze the impact of groundwater drops on land subsidence. Results showed that the highest sinking rate reached 80 mm/year in 2005–2008, and the groundwater decline in the second aquifer influenced the subsidence in the CRFP. Yang, et al. [31] utilized the GPS data to enhance the SBAS results from 2016–2017 in Yunlin County. The results proposed that severe subsidence occurred in Tuku, Yuanchang, and nearby districts, and the surface deformation in the dry season accounted for 60.77% to 73.75% of the total subsidence per year. However, the result only represented one year of observation. Lu, et al. [32] applied the PSI method to analyze SAR images acquired by multiple sensors from 1996–2017 and provide insight into land subsidence issues of the CRFP. The subsiding velocities along the THSR were shortly described in this study, showing the decreasing subsidence rates in recent years. Chen, et al. [33] improved image processing by increasing the number of input SAR images from different satellites over time (1993–2019) and combined the results analyzed by SBAS with other monitoring sensors such as multilayer compaction and groundwater level monitoring wells. Their study provided the historical subsidence in the study area, showing that the largest sinking rates (over 50 mm/year) occurred in the middle-fan section of the CRFP.

While the subsidence patterns in the CRFP have been addressed in several publications, previous studies often simplified InSAR processing by neglecting horizontal movement components, which are critical when assessing the subsidence impact on regional and local scales [33]. In addition, the angular deflections along the THSR derived from InSAR results were occasionally discussed. Therefore, this study aims to extract the recent spatiotemporal development of land subsidence in the CRFP by employing 292 Sentinel-1 SAR images acquired from 2016–2022. The advantages of both PSI and SBAS approaches were applied to produce numerous interferograms by selecting any SAR image pairs satisfying appropriate baseline thresholds. This approach minimizes the decorrelation and topographic error effects, reduces phase aliasing during phase unwrapping [34], and obtains a denser network of measurement points that satisfy amplitude stability and temporal phase coherence conditions. The InSAR-derived line-of-sight (LOS) results were then integrated with CGPS data to extract the vertical displacements, which were subsequently validated by leveling survey data. The distribution of average subsiding velocities in the study area, along with the subsidence profile and angular deflections along the THSR, was presented and discussed based on the calibrated InSAR results. In addition, the study also integrates the InSAR measurements with hydrogeological monitoring data, such as groundwater level, multilayer compaction, and borehole logging records, to propose a straightforward model for quantifying the controlled groundwater level drops between wet and dry seasons. Such allowable groundwater level drops could mitigate the subsidence rate below a specified 40 mm/year goal.

2. Study Area, Datasets, and Methods

2.1. Study Area

The CRFP belongs to the western coastal region of central Taiwan. CRFP has an area of approximately 2500 km², which belongs to Changhua, Yunlin, and the northern part of the Chiayi counties, with a surface elevation ranging from 0 m to 100–150 m above sea level. The CRFP boundary is shaped by Douliu Hill and Bagua Tableland on the eastern side, the Wu River to the north, the Beigang River to the south, and the Taiwan Strait to the west (Figure 1).

The hydrogeological structures of the study area are divided into three sections from the east to the west, including the proximal fan, middle fan, and distal fan [35]. Each section comprised various sedimentary materials, with the average grain sizes decreasing from hilly regions to coastal areas. The sedimentary materials in the CRFP are divided into four primary groups, grading from very coarse grains to very fine grains, which are gravel, coarse sand, fine sand, and clay or silt [35]. These materials are weathered products of rock formations located in the upstream watershed of the study area, such as slate, quartzite, shale, sandstone, and mudstone [12,36]. The borehole profiles suggest that gravel and coarse sand are primarily present at the proximal fan and part of the middle fan, whereas the distal fan mainly witnesses fine-grain materials, including fine sand, clay, and silt (Figure 2).

As it plays an important role in agricultural production areas in Taiwan, the CRFP also demands massive water for irrigation and civilian usage, mainly supplied from groundwater exploitation due to surface water insufficiency in the dry season of a year [13,15,32,37,38]. Based on electricity consumption, Chu, et al. [39] estimated that 13.2 to 26.7 million m³ of groundwater was extracted monthly for irrigation within three townships of Yunlin County (Huwei, Tuku, and Yuanchang) from 2007 to 2017. However, the annual groundwater recharge rates in the alluvial regions are only around 0 to 200 mm/year, much lower than those observed in the eastern hilly areas (800–2000 mm/year) [40,41]. The significant difference between the discharge and recharge of groundwater over a long time was one of the causes resulting in severe land subsidence.

2.2. Geodetic Datasets

2.2.1. Continuous Global Positioning System and Precise Leveling Survey

The CGPS measurements employed in this study were managed by the Institute of Earth Sciences, Academia Sinica (IESAS) [42], with over 400 stations operated by multiple organizations, including IESAS, the Central Weather Bureau (CWB), the Geological Survey and Mining Management Agency, Minstry of Economic Affrairs (GSMMA), and the Ministry of the Interior (MOI). The GNSS data collected in 2006-2023 were processed with GAMIT 10.7/GLOBK 5.35 software packages [43] to obtain GPS position time series. First, double-differenced carrier phase measurements were used to obtain daily solutions with loose constraints. We then used FES2004 (Finite Element Solution tidal mode) [44] for correcting ocean loading and pole-tide. The GMF (Global Mapping Function) [45] was used for the removal of tropospheric delay. The products obtained from GAMIT software, including the variance and covariance matrices of Taiwan and IGS station positions together with other parameters, were subsequently processed with GLOBK to achieve a combined global daily solution in the ITRF2008 reference frame [46].

The precise leveling data were utilized to validate the calibrated InSAR-derived vertical displacements. There are 500 leveling benchmarks in the area of CRFP provided by the Water Resources Agency (WRA). The data measurement is conducted once a year, with the earliest records in the study area dating back to 1992. Leveling surveys were carried out along lines forming loops of hundreds of kilometers. The leveling procedure required that any loop misclosure should be below the criterion of $3\mathrm{m}\mathrm{m}\sqrt{\mathrm{K}}$ (where K is the distance between two adjacent benchmarks in kilometers (km)) [13]. The leveling surveys followed the mean sea level on 1st January 2001 of the Keelung Harbor tide gauge in northern Taiwan [32,47].

2.2.2. SAR Dataset Overview

This study employed 292 SAR images in the ascending track and Interferometric Wide (IW) swath mode [48], acquired by both Sentinel-1A and 1B satellites. The image acquisition time ranged from April 2016 to October 2022. Sentinel-1 is the European Space Agency’s sun-synchronous imaging radar mission (ESA), providing continuous all-weather, day-and-night imagery at C-band radar wavelength. The orbit altitude is around 693 km with an inclination of 98.18°. The satellite couples have the same orbital plane but opposite flight paths and monitor the targeted area together within six days of revisiting. Detailed parameters of the Sentinel-1 data are given in Table 1. The external 30-m Digital Elevation Model (DEM) for topographic phase removal was provided freely by NASA’s Shuttle Radar Topography Mission (SRTM) [49].

2.2.3. SBAS-PSInSAR Processing and Deformation Extraction

Figure 3 shows the SAR image processing workflow for the study. First, pairs of SAR images with perpendicular and temporal baselines less than 200 m and 25 days were coregistered to produce raw interferograms with maximized total correlation. Then, the interferograms were subjected to the Goldstein filter [50] and multilooking by a factor of 3x1 (3 in the range and 1 in the azimuth direction) to lower the noise caused by decorrelation [28]. In addition, the flat-Earth and topographic phases were removed utilizing the 30-m SRTM DEM. Next, persistent scatterers (PSs), or consistently high-coherence pixels, were identified based on amplitude dispersion index (ADI) and phase coherence stability [26,27]. Following that, the Delaunay triangulation method was applied to generate a connected network of all PSs, from which the linear deformation rate and topographic error were extracted through an integration process using the classical growing algorithm for phase unwrapping [51]. The remaining nonlinear deformation rate was separated from the atmospheric perturbation component by applying a low-pass filter in space and a high-pass filter in time [52]. Finally, the atmospheric-free deformation components, including linear and nonlinear deformation rates, from all interferograms were then obtained.

The InSAR-derived results were adjusted utilizing the entire network of GPS stations as multiple reference points [14,32,53]. The process started with projecting the GPS three-dimensional measurements into the Sentinel-1A/B satellites’ line-of-sight (LOS) vector. Next, PSs within a 200 m radius of each GPS station were selected to identify the discrepancies between their average LOS displacements and GPS-derived ones. We assume these discrepancies between GPS stations and the surrounding PSs exist in systematical similarity and could be described by a geostatistical structure. Based on the geostatistical structure, the study employs the Kriging interpolation method to obtain the discrepancy map. All the PS data in the CRFP are then corrected by subtracting the original PS results from the discrepancies. For a specific time of the InSAR result, the repeated process must be conducted by collecting the GPS data for the date of the available SAR image.

2.3. Hydrogeological Datasets

2.3.1. Groundwater Level Dataset Overview

Since 1990, the Water Resources Agency (WRA) has established a dense groundwater level (GWL) monitoring network throughout Taiwan to manage groundwater resources. The network comprises over 350 stations, predominantly located in the western corridor of Taiwan, each station equipped with 1 to 5 wells [54]. These wells provide hourly recordings of groundwater levels from various aquifers at different depths (24 to 307 m). The study area includes 95 GWL monitoring stations, summarized in Table 2. Because the second aquifer in the CRFP is significantly thick, the GWL monitoring wells for this layer are separated into shallow and deep sublayers.

2.3.2. Linear Relation Between Groundwater Level Drops and Subsidence

This study aims to propose a straightforward model to establish thresholds for groundwater extraction to mitigate further land subsidence. The model is based on the correlation between annual groundwater drops and subsidence ranges over the same periods. The yearly groundwater drop is defined as the variance between the highest and lowest groundwater levels recorded annually. Similarly, the subsidence ranges, defined as the differences between the peak and trough values in the corresponding time frames, are derived from the surface deformation data of PSs and leveling benchmarks surrounding each GWL station. Then, linear regression models are applied to fit the datasets in every GWL monitoring well, based on which the expected GWL drop is estimated to keep the subsidence range within a specified threshold of less than 40 mm/year, as proposed by the Water Resources Agency (WRA) of the Ministry of Economic Affairs of Taiwan [55]. Finally, the expected groundwater drops are interpolated to obtain each aquifer’s contour map of thresholds.

2.3.3. Multilayer Compaction Monitoring Wells

Apart from leveling surveys and GPS monitoring for total surface deformation measurements, multilayer compaction monitoring wells (MLCWs) offer detailed insight into the deformation at specific strata, allowing for a comprehensive analysis within the aquifer system. This study employed 32 WRA MLCWs, primarily in Yunlin and partially in Changhua counties. Each well reaches a depth of 300 m and is anchored with 21 to 26 magnetic rings. The positions of these rings were strategically designed based on a combination of hydrogeological models, borehole logging data, and geophysical logging data. Typically, the magnetic rings are installed at critical points, including the boundaries between significant aquifers identified by the Central Geology Survey (CGS), and transitions between fine and coarse sedimentary materials. The depth variations of each magnetic ring are carefully tracked monthly using a probe attached to an indium-alloy measuring tape, ensuring an accuracy of 1 mm. Detailed methodologies for the installation and measurement procedures of the MLCWs can be found in [56].

2.3.4. Estimation of Fine-Grained Sedimentary Material Percentage

Fine-grained sedimentary materials are those in which more than 50% of particles may pass through the No.200 sieve, whose openings are around 74 µm [57]. These materials are further classified as either silt or clay. The percentage of fine-grained sedimentary materials is estimated as the ratio of the aggregate thickness of silt and clay sequences in a borehole to the total drilling length, as follows:

$$P_f={\displaystyle \frac{T_f}{T}}$$

(1)

where $P_f$ stands for the percentage of fine-grained sedimentary materials, $T_f$ is the thickness of fine-grained materials (m), and $T$ is the total drilling length (m).

The borehole logging data used for percentage calculation was from WRA MLCWs and CGS boreholes [58].

3. Results

3.1. SBAS-PSInSAR Average Deformation Map

The spatial distribution of InSAR-derived average deformation rates in the CRFP is illustrated in Figure 4. Over 1.6 million PSs, corresponding to 760 points/km², were identified across the area. Because of their stable backscattering properties, the PSs were mainly detected at structures like houses, buildings, highways, railways, and bridges. Conversely, PSs were scarce in paddy fields, where vegetation or water coverage during rice cropping seasons reduced signal coherence.

The analysis reveals a significant subsidence bowl in Yunlin County, centered in Tuku (TKU) and Yuanzhang (YGZ) districts, with average subsiding velocities between 40 and 60 mm/year, some areas exceeding 70 mm/year. The subsidence rates decrease to around 30 mm/year in the outer regions of this subsidence bowl. Approximately 26% of the study area in Yunlin experienced subsidence rates over 20 mm/year. It is noticed that a 13 km section of the THSR is situated within this severe subsidence zone. In Changhua County, three moderate subsidence bowls in Erlin (ELI), Pitou (PTU), and Xihu (XHU) showed subsidence rates between 20 and 40 mm/year, covering about one-third the size of those in Yunlin.

The InSAR-derived vertical displacements were validated using leveling survey data, the locations of which are shown in Figure 4a. The PSs within a 200m radius of each benchmark were selected, and their average deformation time series was calculated for the validation. Due to the different temporal sampling intervals, validation was performed by comparing the average deformation rates derived from their respective time series. Evaluation metrics, displayed in Figure 4b, include a Mean Absolute Error (MAE) of 4.92 mm/year, a Root-Mean-Square Error (RMSE) of 6.21 mm/year, and a correlation coefficient (r) of 0.91, indicating strong agreement between the measurements from both geodetic techniques.

3.2. The Analysis of Cumulative Compactions

In addition to the InSAR-derived total displacement, which identified areas at risk of severe subsidence, the layer-wise compaction recorded by MLCWs can offer more insights into the underlying mechanism of subsidence. The MLCW measurements are sampled monthly, providing the depths of magnetic rings with respect to the deep foundation of the well (around 300 m depth). The difference in depth recordings at the same ring in different periods is the soil deformation (compaction or expansion) that occurred at the stratigraphic sequence of the ring. The total deformation equals the sum of all deformation values recorded by all rings. Specifically, the sources of groundwater usage might reasonably be inferred based on the depths where major compaction occurred. For example, the presence of major compaction in shallow layers may suggest that the groundwater resource was used for agricultural purposes [56]. On the other hand, if the phenomenon occurs in deep layers, groundwater might reasonably be extracted for domestic and industrial purposes. Figure 5 shows the cumulative compactions recorded by the STES station, which is located near the center of Yunlin’s subsidence bowl. At the first 60-m depth, the primary major compaction can be observed, followed by a secondary major compaction extending to 200 m, corresponding to the highly compressible sand sequences interbedded with fine-grained material layers.

3.3. The Percentage of Fine-grained Materials for the CRFP

An analysis of the cumulative proportion of compaction recorded by MLCWs for every 10-m depth interval revealed that, on average, 50% of cumulative compaction occurred at the first 90-m depth for the entire study area. The proportions increased to 80% and 90% when reaching 170-m and 210-m depths, respectively (see Figure 6a). Based on the analysis, half of the compaction (i.e., 50%) in the study area occurred mainly in the shallow layers at depths less than 100 m. In the CRFP, the available depths of MLCW and geological logging data are inconsistent depending on the practical conditions of the wells. We then used data within the depth of 210 m for the assessment because the data covered more than 90% of the subsidence. Figure 6b shows the contour lines depicting the percentage of fine-grained materials overlapped on the cumulative compactions within 210 m. The percentage of fine-grained sedimentary materials is calculated only at the first 210-m depth of each borehole to demonstrate the distribution of major compaction in the study area. The results show a positive correlation between the high proportion of fine-grained materials and the resulting cumulative compactions. These two main areas are located on the eastern side of Changhua County and in central Yunlin County (Figure 6b). As was recognized by many other cases in other studies, there is high consistency between the cumulative subsidence and fine-grained material proportions in the aquifer systems. The expected high subsidence bowels could commonly show in the areas composed of fine-grained materials. Note that grain sizes less than the silt material type were selected as the fine-grained materials included in the calculations. The clay depositions show a relatively small percentage compared to the silt material based on the borehole loggings obtained from the CRFP.

4. Discussion

4.1. Subsidence Profile Along THSR

The long-term subsidence progress in the CRFP has posed a harmful impact on the THSR section, potentially disrupting operations through technical issues. In addition, land subsidence can weaken the railway-pillar bearing capacity due to negative frictional forces [18,19]. Therefore, for a comprehensive assessment of subsidence impacts on the negative frictional forces that impair THSR operational safety, it is essential to evaluate both vertical and horizontal displacements, as recommended by [33].

Figure 7 shows the cumulative displacements along the THSR section passing through the CRFP study area, as observed at profile AA’ shown in Figure 4. This 49-km railway section is part of a 350-km-long high-speed railway system that runs along the western corridor of Taiwan, from Taipei Station (Taipei City) to Zuoying Station (Kaohsiung City). In Figure 7, the color bar from blue to red visually represents the timeline, with wider gaps between lines indicating higher deformation rates. The subsidence values along the railway were extracted based on the InSAR PS points in the windows of 100-m bandwidth and 50-m distance along the railway. The subsidence profile indicates three remarkably subsiding spots in Xizhou (XZU), Huwei (HWI), and Tuku—Yuangzhang (TKU-YZG) areas, respectively. The sinking spot in Xizhou observed a subsidence rate up to 45 mm/year at the center and around 30 mm/year in the surrounding areas. Towards the south, the subsidence rates significantly increased from HWI to YGZ, with the THSR segments through these areas lying on a massive subsidence bowl, exhibiting sinking rates from 40 to over 60 mm/year. These findings underscore the need to prioritize THSR pillar safety in subsidence assessments and serve as crucial reminders for local officials and researchers. These severe subsidence areas are those shown in Figure 6b. Relatively less subsidence in the XZU than in the HWI could be the relatively small percentage of the fine-grained material.

4.2. The Effects of Angular Deflection on THSR

Land subsidence can lead to the gradual deterioration of buildings and infrastructure over time. The vulnerability of structures to damage caused by subsidence is determined by the degrees of subsidence-related intensity parameters (SRI) commonly utilized in practical geotechnical engineering [59]. An important indicator for assessing potential damage to railways is angular deflection or angular rotation (θ), which can result from differential settlements. Angular deflection at railway pillars is influenced by the compression or rebound of sediment beneath the foundations. The heterogeneity of the material could also influence the vertical variations between pillars. The excessive angular deflection poses a significant risk to the safety of THSR operations. Therefore, the safety code requires that the angular deflection at railway pillars remain under 1/1000 [18,19,60]. The angular deflection at the k^th railway pillar is estimated as follows:

$${\theta }_k=\left|{\displaystyle \frac{d_k-d_{k-1}}{L_n}}+{\displaystyle \frac{d_k-d_{k+1}}{L_m}}\right|$$

(2)

where $d_{k-1}$, $d_k$, and $d_{k+1}$are the displacements at three consecutive pillars: (k-1)^th, k^th, and (k+1)^th. $L_n$ and $L_m$ are the distances from the railway pillar k^th to the (k-1)^th and (k+1)^th pillars, respectively.

Traditionally, angular deflections along the THSR are estimated based on yearly leveling survey measurements with low measurement density and limited temporal resolution. In the study, utilizing the PSs obtained from the SBAS-PSInSAR process can capture surface deformations in high temporal resolution, which enables the calculation of pseudo-real-time values for angular deflections, providing a more accurate and detailed assessment of the risk and potential damage through relatively short periods [61]. This study considered that the inter-pillar distances are approximately 50 m for the analysis.

The variations in angular deflections over time are derived from high-temporal measurements of differential settlement, and results could provide an assessment of potential damage to railway pillars. This approach differs from typical leveling survey calculations, which analyze angular deflections at a specific time. Figure 8 represents the box plot analysis of angular deflections calculated from the cumulative displacements at the end of the research period and categorized based on the districts crossed by the railway system. The analysis shows that the angular deflections in Shetou, Xizhou, Tuku, and Yuanzhang districts are higher than those in the other areas, with plenty of pillars approaching the safety code. This observation suggests that the THSR segments in these places are more likely to suffer damage induced by nonuniform vertical variations (i.e., the subsidence or rebound). In addition to the spatial viewpoint, a temporal perspective of angular deflections can be achieved by analyzing the angular deflections calculated from differential displacements between successive SAR image acquisition dates. Figure 9 shows that the values exceeding the 90th percentile of angular deflections along the railway occurred primarily between April and October. The result indicates that high angular deflections are commonly observed in the wet seasons, which are not the period with significant subsidence in dry seasons. It is noted that the formula of angular deflection in Equation (1) implies that this parameter demonstrates the extent of non-uniformity in surface deformation rather than subsidence magnitude. The possible factors that control the nonuniform deformation might be the heterogeneity of the material under the foundations of the railway. This heterogeneity could be the natural deposition of the aquifer materials and the soil reinforcement for the railway construction. In summary, the spatiotemporal analysis of angular deflections derived from PSs shed light on the THSR segments susceptible to nonuniform surface deformation and when the non-uniformity occurred at its peak. For safety issues, more angular deflection measurements in the wet seasons might be required.

4.3. Proposing Groundwater Drop Threshold for Subsidence Mitigation

A large number of GWL monitoring stations and leveling benchmarks have been installed in the CRFP. The study integrated the groundwater monitoring data with the PSs identified during InSAR processing. The data integration was implemented by constructing linear regression models between annual groundwater drops and subsidence magnitudes, as mentioned in Section 3.2. Due to seasonal effects on the variations of groundwater levels and subsidence [62], the time series of GWL monitoring and surface deformations exhibit periodic patterns. Therefore, the annual groundwater drops and the subsidence amplitudes are equal to the peak-to-peak values of the sinusoidal waveform. In the study, the InSAR provides maximal subsidence between wet and dry seasons, marked as total subsidence near the groundwater monitoring wells. The groundwater level drops corresponding to the total subsidence are identified based on the multilevel groundwater monitoring network in different aquifers in the CRFP. For an observation well at a specific aquifer depth, the relationship between total subsidence and groundwater level drops is derived based on the data from 2016 to 2023. With a specified subsidence mitigation target, for example, 30 or 40 mm/year at the monitoring well location, the allowable GWL drops could be obtained based on the relationship. Figure 10 shows the estimated GWL drops for each aquifer to mitigate subsidence amplitudes under 40 mm/year, overlapped by the borders of areas suffering from severe subsidence. The warmer colors indicate lower thresholds of GWL drops and imply that the annual GWL drops in a specific layer are significantly sensitive to the total subsidence. With this perspective, aquifers 1, 3, and 4 from shallow to deep may need more attention for subsidence mitigation, and the controls of the GWL drops in these layers might be critical. Note that aquifer 2 in the study is further divided into two sub-aquifers named aquifer 2-shallow and aquifer 2-deep.

In the first aquifer layer, the GWL drop thresholds, ranging from 3 to 5 m, correspond to the areas prone to severe subsidence and their vicinity, implying substantial GWL drops in this aquifer remarkably impact the total subsidence. The first aquifer is the shallowest layer and is known as the cost-effective source of irrigation for agricultural activities. Note that we had specified the subsidence rate of 40 mm/year for the study area. Other desired subsidence rates could follow the same strategies to quantify the controlled GWL drops for different aquifers. A similar range of GWL drop thresholds can be observed in the deepest aquifer, the source of industrial and domestic water usage. However, the estimated GWL drops in Aquifer 4 may be uncertain due to the limited GWL monitoring wells. Aquifer 3 is also worth consideration due to the moderate GWL drop thresholds, which range from 4 to 6 m and primarily coincide with the subsidence-prone areas. The GWL drop thresholds observed in aquifers 3 and 4 remind the local authorities to pay attention to the amount of groundwater usage for domestic and industrial purposes. On the other hand, urgent actions are required to manage the groundwater usage of the shallowest aquifer since the estimated GWL drops for this layer are remarkably low, and half of the major compaction in the study area occurred at the shallow sequences, as shown in Figure 6.

5. Conclusions

This study employed geodetical and hydrogeological data sets to assess the land subsidence evolution in the CRFP and suggest a straightforward model for subsidence mitigation. The data analysis of WRA MLCWs and CGS boreholes suggested that 50% of major compaction occurred within the first 90-m depth, corresponding to the shallow layers of the study area. Additionally, the major compaction positively correlated with the fine-grained material proportion. The InSAR results showed a prominent subsidence bowl centered in Yunlin County, with average subsidence velocities around 40 to 60 mm/year. Three medium-sized subsidence bowls were also detected in Changhua County, with subsiding rates varying 20 to 40 mm/year. The InSAR results agreed well with precise leveling data, as proven by three evaluation metrics (MAE = 4.92 mm/year; RMSE = 6.21 mm/year; and r = 0.91). The subsidence profile along the THSR revealed multiple sinking spots in the XZU, HWI, TKU, and YZG districts. In addition, the analysis of angular deflections identified railway segments highly prone to nonuniform subsiding rates. These findings could bring the attention of the local officials to the subsidence escalation in these areas and urge them to introduce appropriate subsidence-control policies. The model developed based on the linear relationship between the annual records of groundwater level drops and subsidence amplitudes suggested that aquifers 1, 3, and 4 require more attention in terms of subsidence mitigation. Specifically, to maintain the subsidence rate under 40 mm/year, the GWL drop should be limited to 3–5 m for the shallowest aquifer and below 4–6 m for aquifers 3 and 4. This study demonstrated the reliability of fusing the InSAR results with traditional CGPS and leveling survey data and the potential of integrating geodetic and hydrogeological monitoring data for further subsidence mitigation research.