Land Subsidence-Induced Horizontal Displacement Along the High-Speed Rail in Central Taiwan: An Integrated Multi-Temporal InSAR, GNSS, and Leveling Approach

1. Introduction

Land subsidence driven by excessive groundwater extraction is a pervasive geohazard in coastal alluvial plains worldwide [1,2]. In the Choushui River alluvial fan of central Taiwan (Figure 1), decades of intensive groundwater pumping for agricultural irrigation and aquaculture have resulted in significant surface subsidence, with cumulative settlement exceeding 2 m in some areas since the 1970s [3,4]. This ongoing deformation poses a critical threat to the Taiwan High-Speed Rail (THSR), which traverses the subsidence-affected Changhua and Yunlin counties over a distance of approximately 40 km.

Subsidence monitoring has traditionally relied on precise leveling surveys and continuous GPS (CGPS) networks, which provide point-wise measurements with millimeter-level accuracy [5,6]. However, these techniques are inherently limited in spatial coverage. Interferometric Synthetic Aperture Radar (InSAR), particularly multi-temporal InSAR (MTI) methods such as Persistent Scatterer InSAR (PSI) and Small Baseline Subset (SBAS), offer spatially dense deformation measurements with sub-centimeter precision over large areas [7,8]. InSAR has been widely applied to subsidence monitoring in alluvial settings globally [9,10].

Previous studies have documented the subsidence patterns in the Choushui River alluvial fan using various geodetic techniques. Liu et al. [3] characterized land subsidence using geological and hydrogeological data. Chang et al. [11] investigated the relationship between seasonal groundwater fluctuations and surface deformation using SAR interferometry. Hung et al. [12] employed multiple sensors including InSAR, GPS, and leveling for monitoring severe aquifer-system compaction in the southern fan. Lu et al. [13] examined the relationship between surface displacement and groundwater level changes using remote sensing approaches.

While most studies focus on the vertical component of subsidence, theoretical and observational evidence indicates that aquifer compaction also induces significant horizontal displacements directed toward subsidence centers [14,15]. Burbey [14] demonstrated that land subsidence is accompanied by horizontal displacement in basin-fill deposits through the Poisson effect during vertical compaction. Bawden et al. [16] used InSAR and GPS to detect horizontal motions associated with groundwater pumping in the Los Angeles Basin. These horizontal motions are critical for infrastructure integrity, particularly for linear structures such as railways, where differential lateral displacements can cause track misalignment and structural distress.

Isolating horizontal motions from InSAR remains challenging because InSAR is primarily sensitive to the line-of-sight (LOS) direction. The 2.5D decomposition of ascending and descending LOS velocities yields east-west and vertical components, but the north-south component is inherently unresolvable due to the near-polar orbit geometry of SAR satellites. Furthermore, the decomposed EW and vertical fields contain systematic biases from atmospheric phase delays, orbital ramps, and imperfect decomposition geometry that must be corrected using independent ground-truth data.

The primary contributions of this study are: (1) the development of an integrated GNSS and leveling calibration framework for InSAR-derived EW and vertical velocity fields that dramatically improves accuracy; (2) the comprehensive characterization of both vertical and horizontal deformation fields along the THSR corridor with unprecedented spatial density; (3) the demonstration that GNSS-derived north-south velocities provide critical complementary information in the InSAR blind direction; and (4) a quantitative assessment of the potential hazards posed by subsidence-induced horizontal displacement to the THSR structure.

2. Study Area

The study area encompasses Changhua and Yunlin counties in central western Taiwan, situated on the Choushui River alluvial fan, one of the largest alluvial fans in Taiwan (Figure 1, Figure 2). The fan extends approximately 50 km from the Western Foothills to the Taiwan Strait coastline, covering an area of roughly 2,000 km².

The Quaternary alluvial deposits form a complex multi-layered aquifer system consisting of alternating gravel, sand, and clay layers extending to depths exceeding 300 m (Figure 3). The hydrostratigraphy is characterized by four principal aquifer units (Aquifer I through IV) separated by aquitard layers of marine clay and silt [3,4]. Previous studies indicate that most groundwater resources come from the second and third aquifers.

Groundwater extraction intensified from the 1970s onward, primarily for agriculture, aquaculture, and industrial use. The consequent decline in piezometric head has caused compaction of the compressible clay aquitard layers through effective stress increase. This compaction is predominantly inelastic when effective stress exceeds the preconsolidation stress, resulting in permanent subsidence. The rate of anthropogenic subsidence far exceeds the natural subsidence rate of 0.6–2.1 mm/yr estimated from borehole data [17].

In Changhua County, three main subsidence centers are located near Xihu, Erlin, and Xizhou. In Yunlin County, the subsidence centers are near Baozhong, Tuku, Huwei, and Yuanchang [13]. The maximum subsidence rates in Tuku reach 60–67 mm/yr. The THSR traverses both the Xizhou (Changhua) and Tuku (Yunlin) subsidence zones, making it particularly vulnerable to both vertical and horizontal deformation.

3. Materials and Methods

The overall research methodology and workflow adopted in this study is illustrated in Figure 4. The processing chain consists of four main stages: (1) SAR data acquisition and InSAR time-series processing using the SBAS/MTI approach; (2) 2.5D decomposition of ascending and descending LOS velocities into east-west and vertical components; (3) calibration of the decomposed velocity fields using multi-source ground-truth data including GNSS and precise leveling observations through Ordinary Kriging interpolation and polynomial trend surface fitting; and (4) integrated analysis of the calibrated three-dimensional velocity field to assess subsidence-induced horizontal deformation and its implications for the THSR corridor. The following subsections describe each data source and processing step in detail.

3.1. SAR Data and 2.5D Decomposition

Sentinel-1A/B C-band SAR data were acquired from ascending track A69 and descending track D105 over the period 2015–2021 (Table 1). SAR data were processed using the SBAS/MTI approach implemented in MintPy, which exploits small spatial and temporal baselines to minimize decorrelation while maximizing temporal sampling density. The processing chain included co-registration, interferogram generation, topographic phase removal using the SRTM 1-arc-second DEM, phase unwrapping using the minimum cost flow algorithm, and time-series inversion.

The ascending and descending LOS velocity fields were decomposed into east-west (EW) and vertical components using 2.5D decomposition [18]. Given the near-polar orbit geometry of Sentinel-1, InSAR is inherently insensitive to the north-south component, which is neglected in the decomposition (V_N ≈ 0). The system of equations for each co-located pixel is:

$$V_{LOS}^{asc}=-sin\left({\theta }_{asc}\right)cos\left({\alpha }_{asc}\right)V_{EW}+cos\left({\theta }_{asc}\right)V_U$$

(1)

$$V_{LOS}^{desc}=-sin\left({\theta }_{desc}\right)cos\left({\alpha }_{desc}\right)V_{EW}+cos\left({\theta }_{desc}\right)V_U$$

(2)

where θ is the radar incidence angle and α is the satellite heading angle. The resulting EW and vertical velocity fields serve as the starting point for the subsequent GNSS and leveling calibrations described in this study.

3.2. GNSS Data

Two categories of GNSS observations were incorporated. First, 44 continuous GPS (CGPS) stations from the Taiwan CGPS network operated by Academia Sinica and the Central Geological Survey provide daily position solutions in the ITRF2014 reference frame. Second, 42 campaign GNSS stations surveyed by the Water Resources Agency and the Industrial Technology Research Institute densify the spatial coverage within the subsidence zone. All GNSS velocities are expressed relative to station GS26 (120.6460° E, 23.8698° N), located on stable ground east of the subsidence zone.

The GNSS velocity field provides three-component (east, north, up) ground-truth for calibrating the InSAR-derived EW and vertical fields. Importantly, the north-south component from GNSS provides information in the direction that InSAR inherently cannot detect.

3.3. Leveling Data

Precise leveling data were obtained from the Water Resources Agency, encompassing annual survey campaigns from 2016 to 2021 (ROC years 105–110). A total of 922 leveling benchmarks were used, covering Changhua (368 benchmarks) and Yunlin (554 benchmarks) counties. The leveling surveys were conducted each May, providing six elevation snapshots from which vertical velocities were derived by linear regression. Yunlin benchmarks exhibit a median subsidence rate of −17.9 mm/yr with maximum exceeding −72.5 mm/yr, while Changhua shows a median of −7.1 mm/yr. The leveling observations provide sub-millimeter accuracy (Class II leveling, ~0.5–1 mm/yr) and serve as the primary reference for calibrating the InSAR vertical displacement field.

3.4. East-West Velocity Field Correction with GNSS Kriging

The InSAR-derived EW velocity field from the 2.5D decomposition contains spatially correlated systematic biases arising from atmospheric phase residuals, orbital ramps, and imperfect decomposition geometry. These errors manifest as poor correlation with GNSS EW velocities: initial comparison yields R² = 0.147 with RMSE = 4.24 mm/yr (n = 89 stations) (Figure 5, Figure 8a). At the continuous GNSS stations, the correlation is R² = 0.207 (RMSE = 4.66 mm/yr, n = 47), while the campaign GNSS stations show essentially no correlation (R² = 0.003, n = 42), indicating that the systematic biases dominate the true EW deformation signal. The mean bias (InSAR minus GNSS) before correction is −1.31 mm/yr for all stations combined, with the continuous GNSS stations showing a near-zero bias of −0.18 mm/yr and the campaign GNSS stations exhibiting a larger negative bias of −2.57 mm/yr. After Kriging correction using all GNSS stations, the bias is virtually eliminated to +0.01 mm/yr (continuous: +0.00 mm/yr; campaign: +0.02 mm/yr), confirming that the Kriging interpolation effectively removes both the spatially correlated systematic errors and the mean offset of the EW velocity field (Table 2).

To correct these biases, we computed the EW velocity residuals at each GNSS station:

$$\Delta V_{EW,i}=V_{EW,i}^{GNSS}-V_{EW,i}^{InSAR}$$

(3)

These residuals were spatially interpolated over the entire study area using Ordinary Kriging with a spherical variogram model. The Kriging-interpolated correction surface was then added to the original InSAR EW velocity field. The correction was applied using both CGPS and campaign GNSS stations combined. After correction, the agreement improves dramatically to R² = 0.992 with RMSE = 0.24 mm/yr (Figure 7, Figure 8b). The improvement demonstrates the critical importance of dense ground-truth observations for calibrating InSAR-derived horizontal velocity fields (Figure 7).

Figure 6. GNSS Kriging correction field (difference surface).

Figure 6. GNSS Kriging correction field (difference surface).

Figure 7. Corrected InSAR EW velocity field after GNSS Kriging calibration.

Figure 7. Corrected InSAR EW velocity field after GNSS Kriging calibration.

3.5. Vertical Velocity Field Correction with Leveling and GNSS

The InSAR-derived vertical (UD) velocity field from the 2.5D decomposition similarly contains long-wavelength biases (Figure 11). Unlike the EW component where Kriging interpolation of point residuals is applied, the vertical field is corrected using a 2nd-degree polynomial trend surface fitted to the combined leveling (Figure 9) and GNSS control points. This approach is more appropriate for the vertical component because: (1) the 922 leveling benchmarks provide sufficient redundancy to define the spatial trend; (2) the residuals primarily reflect smooth, low-frequency atmospheric and orbital biases; and (3) applying Kriging to sparse GNSS data risks over-fitting local noise.

Before correction, the InSAR vertical velocity field exhibits varying degrees of agreement depending on the validation dataset (Figure 13, left; Table 2): the continuous GNSS stations show RMSE = 11.56 mm/yr (R² = 0.754) with a systematic positive bias of +5.05 mm/yr (n = 44), the campaign GNSS stations show RMSE = 5.60 mm/yr (R² = 0.820) with a negative bias of −3.67 mm/yr (n = 42), and the leveling benchmarks show RMSE = 6.03 mm/yr (R² = 0.889) with a positive bias of +3.47 mm/yr (n = 788). The contrasting bias signs between continuous GNSS (+5.05 mm/yr) and campaign GNSS (−3.67 mm/yr) suggest that the two datasets are affected differently by reference frame alignment and temporal coverage differences.

After the polynomial trend surface correction (Figure 12), the leveling validation improves to RMSE = 4.85 mm/yr (R² = 0.898, n = 907) with the bias virtually eliminated to −0.47 mm/yr, representing a 17% RMSE reduction and an 86% bias reduction (Table 2). The continuous GNSS bias decreases from +5.05 to +1.53 mm/yr (70% reduction). However, the campaign GNSS shows degraded performance after correction (RMSE from 5.60 to 7.70 mm/yr, bias from −3.67 to −6.34 mm/yr) (Figure 13, right; Table 2), which is attributed to short occupation periods (1–3 days) and incomplete temporal overlap with the InSAR observation window. Given this limitation, the leveling-based statistics (n = 907, RMSE = 4.85 mm/yr, R² = 0.898) are considered the most reliable external accuracy measure for the corrected vertical velocity field (Table 2).

The polynomial trend surface is fitted using weighted least squares:

$$T\left(x,y\right)={\beta }_0+{\beta }_{1}x+{\beta }_{2}y+{\beta }_3{x}^2+{\beta }_{4}xy+{\beta }_5{y}^2$$

(4)

where x and y are normalized geographic coordinates, and the weights are assigned as w = 3 for GNSS stations (continuous observations with consistent reference frame) and w = 1 for leveling benchmarks (annual discrete observations). The total control point count is 918 leveling benchmarks + 38 GNSS stations = 956 points. The fitted trend surface ranges from −18.84 to +3.97 mm/yr and is subtracted from the original UD field to produce the corrected vertical velocity field.

Table 2 summarizes the calibration results for both EW and vertical components. Bias is defined as the mean difference (InSAR minus ground-truth), representing the systematic offset of the velocity field. The EW correction achieves near-perfect agreement (R² = 0.992) because the same GNSS observations constrain the Kriging field. For the vertical component, the bias provides critical additional information beyond RMSE and R²: the leveling bias decreases from +3.47 to −0.47 mm/yr, indicating that the polynomial trend surface virtually eliminates the systematic offset. Similarly, the continuous GNSS bias decreases from +5.05 to +1.53 mm/yr (a 70% reduction). However, the campaign GNSS bias worsens from −3.67 to −6.34 mm/yr, likely due to short occupation times (1–3 days) and incomplete temporal overlap with the InSAR observation period. The leveling-based RMSE of 4.85 mm/yr (n = 907) represents the most reliable external accuracy estimate for the corrected vertical velocity field.

4. Results

4.1. Calibrated East-West Velocity Field

The InSAR-derived EW velocity field before and after GNSS Kriging correction is presented in Figure 8. Before correction, the EW field contains a pronounced long-wavelength artifact that obscures the true deformation signal, with values ranging from approximately −8 to +8 mm/yr in a spatially incoherent pattern. The Kriging correction field (Figure 8b) captures this systematic bias, which varies smoothly across the study area with amplitudes of several mm/yr.

After correction (Figure 8c), the EW velocity field reveals a clear and physically meaningful pattern associated with subsidence-induced horizontal deformation. In the subsidence zone, the corrected EW velocities show convergent motion directed toward the subsidence centers: eastward motion on the western side and westward motion on the eastern side. The maximum corrected EW velocities reach approximately 5–8 mm/yr near the steepest gradient of the subsidence bowl. The scatter plot comparison confirms the dramatic improvement from R² = 0.147 to R² = 0.992 (Figure 8d).

Figure 8. Scatter plots comparing InSAR EW velocity with GNSS EW velocity under different correction scenarios. Blue circles represent continuous GNSS stations (n = 47); red triangles represent campaign GNSS stations (n = 42). The dashed line is the 1:1 reference; the solid line is the linear regression fit. (a) Before correction: all 89 stations show poor agreement (R² = 0.147, RMSE = 4.24 mm/yr, Bias = −1.31 mm/yr), with large scatter and a regression slope of 0.62, indicating that the InSAR EW field is dominated by systematic biases. (b) After Kriging correction using continuous GNSS stations only: the overall fit improves substantially (R² = 0.815, RMSE = 1.15 mm/yr, Bias = −0.18 mm/yr), but the campaign GNSS stations (red triangles) still deviate from the 1:1 line, suggesting that the 47 continuous stations alone do not provide sufficient spatial density to fully resolve the correction field. (c) Same as (a), shown for comparison with (d). (d) After Kriging correction using both continuous and campaign GNSS stations: near-perfect agreement is achieved (R² = 0.992, RMSE = 0.24 mm/yr, Bias = +0.01 mm/yr), with both continuous and campaign stations tightly clustered along the 1:1 line and a regression slope of 0.97. The comparison between (b) and (d) demonstrates that incorporating the 42 campaign GNSS stations is critical for achieving high-accuracy EW velocity field calibration, as they fill the spatial gaps in the continuous network within the subsidence zone.

Figure 8. Scatter plots comparing InSAR EW velocity with GNSS EW velocity under different correction scenarios. Blue circles represent continuous GNSS stations (n = 47); red triangles represent campaign GNSS stations (n = 42). The dashed line is the 1:1 reference; the solid line is the linear regression fit. (a) Before correction: all 89 stations show poor agreement (R² = 0.147, RMSE = 4.24 mm/yr, Bias = −1.31 mm/yr), with large scatter and a regression slope of 0.62, indicating that the InSAR EW field is dominated by systematic biases. (b) After Kriging correction using continuous GNSS stations only: the overall fit improves substantially (R² = 0.815, RMSE = 1.15 mm/yr, Bias = −0.18 mm/yr), but the campaign GNSS stations (red triangles) still deviate from the 1:1 line, suggesting that the 47 continuous stations alone do not provide sufficient spatial density to fully resolve the correction field. (c) Same as (a), shown for comparison with (d). (d) After Kriging correction using both continuous and campaign GNSS stations: near-perfect agreement is achieved (R² = 0.992, RMSE = 0.24 mm/yr, Bias = +0.01 mm/yr), with both continuous and campaign stations tightly clustered along the 1:1 line and a regression slope of 0.97. The comparison between (b) and (d) demonstrates that incorporating the 42 campaign GNSS stations is critical for achieving high-accuracy EW velocity field calibration, as they fill the spatial gaps in the continuous network within the subsidence zone.

4.2. Calibrated Vertical Velocity Field

The vertical velocity field calibration is shown in Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13. The leveling benchmark velocity distribution (Figure 9) provides the spatial reference for the correction. Before correction (Figure 10), the InSAR vertical field captures the overall subsidence pattern but exhibits a systematic positive bias of +3.47 mm/yr relative to the leveling data, arising from atmospheric artifacts and orbital residuals propagated through the 2.5D decomposition.

The polynomial trend surface correction field (Figure 11) ranges from −18.84 to +3.97 mm/yr, indicating that the long-wavelength bias varies significantly across the study area. After correction (Figure 12), the vertical field shows improved agreement with ground-truth data. The maximum subsidence rate in the Yunlin Tuku area exceeds 60 mm/yr, and three secondary subsidence centers in Changhua exhibit rates of 15–35 mm/yr. The subsidence signal diminishes rapidly toward the piedmont zone in the east. The quantitative improvement is evaluated by comparing InSAR vertical velocities against three independent datasets before and after correction (Figure 13): the leveling benchmarks show RMSE improvement from 5.98 to 5.09 mm/yr (R² from 0.879 to 0.889) with the systematic bias reduced from +3.27 to −0.66 mm/yr (n = 918); the continuous GNSS bias decreases from +5.05 to +1.53 mm/yr (RMSE = 11.56 to 10.27 mm/yr, n = 44); however, the campaign GNSS RMSE degrades from 5.60 to 7.70 mm/yr (n = 42), attributed to short occupation periods and temporal mismatch with the InSAR observation window.

The leveling time series analysis (2016–2021) confirms that subsidence was temporally stable during the study period, with median rates of −17.9 mm/yr for Yunlin (551 benchmarks) and −7.1 mm/yr for Changhua (368 benchmarks). The near-linear temporal trends at most benchmarks justify the comparison of InSAR time-averaged rates with leveling annual-average velocities.

Figure 9. Vertical velocity field calibration. (a) Leveling benchmark velocity distribution across Changhua and Yunlin counties.

Figure 9. Vertical velocity field calibration. (a) Leveling benchmark velocity distribution across Changhua and Yunlin counties.

Figure 10. InSAR vertical velocity before polynomial trend surface correction.

Figure 10. InSAR vertical velocity before polynomial trend surface correction.

Figure 11. Polynomial trend surface correction field (mm/yr).

Figure 11. Polynomial trend surface correction field (mm/yr).

Figure 12. Corrected InSAR vertical velocity field.

Figure 12. Corrected InSAR vertical velocity field.

Figure 13. Scatter plots comparing InSAR vertical (UD) velocity with ground-truth measurements before (left) and after (right) polynomial trend surface correction. Three categories of validation data are shown: blue circles for continuous GNSS stations (Cont. GNSS), red triangles for campaign GNSS stations (Camp. GNSS), and green squares for leveling benchmarks (Leveling). The dashed line is the 1:1 reference; colored solid lines are the linear regression fits for each dataset. (Left) Before correction (InSAR 2.5D, GNSS-corrected only): the continuous GNSS stations show RMSE = 11.56 mm/yr (R² = 0.754, Bias = +5.05 mm/yr, n = 44), with a notable systematic positive offset indicating that the InSAR UD field overestimates the vertical velocity relative to GNSS; the campaign GNSS stations show better agreement with RMSE = 5.60 mm/yr (R² = 0.820, Bias = −3.67 mm/yr, n = 35), though with a negative bias of opposite sign; the leveling benchmarks provide the densest comparison with RMSE = 5.98 mm/yr (R² = 0.879, Bias = +3.27 mm/yr, n = 918), showing good overall correlation but a clear positive offset from the 1:1 line. (Right) After correction (GNSS + Leveling polynomial trend surface): the continuous GNSS bias is substantially reduced from +5.05 to +1.53 mm/yr (RMSE = 10.27 mm/yr, R² = 0.769, n = 44); the leveling comparison improves to RMSE = 5.09 mm/yr (R² = 0.889, Bias = −0.66 mm/yr, n = 918), with the data points tightening around the 1:1 line; however, the campaign GNSS performance degrades after correction (RMSE from 5.60 to 7.70 mm/yr, Bias from −3.67 to −6.34 mm/yr, n = 42), attributed to short occupation periods and temporal mismatch with the InSAR observation window. The contrasting responses of the three datasets to the correction highlight the importance of using temporally consistent, spatially dense observations (i.e., leveling benchmarks) as the primary validation reference for vertical velocity field calibration.

Figure 13. Scatter plots comparing InSAR vertical (UD) velocity with ground-truth measurements before (left) and after (right) polynomial trend surface correction. Three categories of validation data are shown: blue circles for continuous GNSS stations (Cont. GNSS), red triangles for campaign GNSS stations (Camp. GNSS), and green squares for leveling benchmarks (Leveling). The dashed line is the 1:1 reference; colored solid lines are the linear regression fits for each dataset. (Left) Before correction (InSAR 2.5D, GNSS-corrected only): the continuous GNSS stations show RMSE = 11.56 mm/yr (R² = 0.754, Bias = +5.05 mm/yr, n = 44), with a notable systematic positive offset indicating that the InSAR UD field overestimates the vertical velocity relative to GNSS; the campaign GNSS stations show better agreement with RMSE = 5.60 mm/yr (R² = 0.820, Bias = −3.67 mm/yr, n = 35), though with a negative bias of opposite sign; the leveling benchmarks provide the densest comparison with RMSE = 5.98 mm/yr (R² = 0.879, Bias = +3.27 mm/yr, n = 918), showing good overall correlation but a clear positive offset from the 1:1 line. (Right) After correction (GNSS + Leveling polynomial trend surface): the continuous GNSS bias is substantially reduced from +5.05 to +1.53 mm/yr (RMSE = 10.27 mm/yr, R² = 0.769, n = 44); the leveling comparison improves to RMSE = 5.09 mm/yr (R² = 0.889, Bias = −0.66 mm/yr, n = 918), with the data points tightening around the 1:1 line; however, the campaign GNSS performance degrades after correction (RMSE from 5.60 to 7.70 mm/yr, Bias from −3.67 to −6.34 mm/yr, n = 42), attributed to short occupation periods and temporal mismatch with the InSAR observation window. The contrasting responses of the three datasets to the correction highlight the importance of using temporally consistent, spatially dense observations (i.e., leveling benchmarks) as the primary validation reference for vertical velocity field calibration.

4.3. Horizontal Velocity Field and Subsidence Centers

The horizontal velocity field, derived from the calibrated InSAR EW component and supplementary GNSS north-south observations, is superimposed on the vertical subsidence rate in Figure 7. The background color map displays the corrected InSAR vertical velocity field (color scale: −45 to +5 mm/yr), while the arrows represent the combined horizontal velocity vectors referenced to the Yunlin subsidence center (red star, ~120.24° E, 23.36° N). The figure reveals a prominent convergent horizontal velocity pattern centered on the major subsidence zones: arrows from all directions point radially inward toward the subsidence center, with magnitudes of 2–10 mm/yr. The maximum horizontal velocities occur not at the subsidence center itself, but along the steepest gradient of the vertical subsidence bowl, approximately 5–15 km from the center, where the transition from rapid subsidence to stable ground generates the strongest lateral strain. In the Tuku area (Yunlin), where vertical subsidence exceeds 60 mm/yr, the surrounding horizontal velocities reach 8–10 mm/yr; in the Xizhou area (Changhua), with vertical subsidence of 15–35 mm/yr, horizontal velocities of 3–5 mm/yr are observed. The spatial correspondence between the subsidence magnitude and the horizontal velocity amplitude provides direct observational evidence that horizontal displacement is mechanically coupled to vertical compaction.

This convergent pattern is consistent with the Poisson effect associated with vertical compaction of the aquifer system [14,16]. As the compressible clay aquitard layers undergo vertical compression due to effective stress increase, the lateral contraction induced by the Poisson ratio (ν ≈ 0.3–0.4 for clay) generates horizontal strain directed toward the center of compression. This strain is transmitted upward through the overlying soil column to the surface, producing the observed radially convergent motion. The spatial gradient of horizontal velocity provides a direct measure of the horizontal strain rate. In the vicinity of the Tuku subsidence center, where the horizontal velocity changes from approximately +5 mm/yr (eastward) to −5 mm/yr (westward) over a distance of roughly 10 km, the resulting strain rate reaches approximately 10⁻⁶/yr. This strain rate, integrated over the multi-decadal subsidence history, implies cumulative horizontal shortening on the order of several centimeters across the subsidence zone. Furthermore, the asymmetry of the arrow pattern with stronger convergence from the western (coastal) side than from the eastern (piedmont) side reflects the asymmetric geometry of the subsidence bowl, which extends more broadly toward the coast due to the thicker compressible sediments in the distal-fan zone.

For an almost north-south oriented structure such as the THSR, the east-west component of horizontal deformation is particularly hazardous because it acts perpendicular to the rail alignment. As shown in Figure 14, the THSR corridor (running approximately N-S through the study area) traverses the transition zone where the horizontal velocity arrows rotate from eastward on the northern flank to westward on the southern flank of the subsidence bowl, creating a differential EW displacement across the rail alignment. This differential displacement can cause: (1) lateral misalignment of the track, compromising ride quality and safety at operating speeds of up to 300 km/h; (2) additional bending stresses on viaduct piers and foundations, as adjacent piers experience different magnitudes and directions of horizontal forcing; (3) cumulative damage to expansion joints and bearing systems, which are designed to accommodate thermal expansion but not sustained lateral drift; and (4) progressive tilting of viaduct columns in the subsidence transition zone, where combined vertical settlement and horizontal displacement impose complex three-dimensional loading on the foundation system.

4.4. THSR Corridor Profile Analysis

Velocity profiles extracted along the THSR alignment reveal two principal zones of concern (Figure 15). The Xizhou zone in southern Changhua (distance ~36 km from the northern end of the study area, lat ~23.84° N) shows vertical subsidence rates of 15–25 mm/yr with an EW velocity gradient of approximately 3–5 mm/yr over a 5-km span. The Tuku zone in Yunlin (distance ~62 km, lat ~23.65° N) exhibits more severe conditions, with vertical subsidence exceeding 50 mm/yr and EW velocity change of up to 8 mm/yr across the subsidence transition zone.

The EW velocity profile along the THSR shows a characteristic pattern across each subsidence center: the velocity increases as the profile approaches the center, then undergoes a transition across it, reflecting the convergent horizontal strain field. The maximum differential EW velocity along the THSR reaches approximately 10 mm/yr over a baseline of 10 km, corresponding to a horizontal strain rate of approximately 10⁻⁶/yr. The combination of vertical differential settlement and lateral displacement creates a three-dimensional strain field that significantly complicates structural health monitoring and maintenance planning for the THSR.

Figure 16 reveals a clear anti-correlation between the vertical subsidence rate and the east-west horizontal velocity along the THSR alignment. At the Tuku subsidence center (distance ~62 km), the vertical velocity reaches its maximum (~60 mm/yr downward), while the EW velocity crosses zero, transitioning from eastward motion (positive) on the northern flank to westward motion (negative) on the southern flank. This zero-crossing at the subsidence center, flanked by opposing EW velocities, is the hallmark of convergent horizontal strain driven by vertical compaction [14].

At the Xizhou subsidence center (distance ~36 km), a similar but weaker pattern is observed: vertical subsidence of 15–25 mm/yr is accompanied by a localized EW velocity gradient of approximately 5 mm/yr. Between the two subsidence centers (distance 36–62 km), the EW velocity maintains a relatively stable eastward value of 5–10 mm/yr, indicating that this inter-center segment experiences translational rather than differential horizontal motion.

The spatial correlation between the subsidence bowl geometry and the horizontal velocity pattern has important engineering implications. The maximum differential horizontal displacement does not occur at the subsidence center itself, but rather at the inflection points on the flanks of the subsidence bowl, where the vertical velocity gradient is steepest. For the THSR, this means that the track segments most vulnerable to lateral misalignment are located 5–10 km north and south of each subsidence center, rather than directly above it. The cumulative effect over the 50-year design life translates to differential horizontal displacements of 15–40 cm at these critical transition zones, which may exceed the design tolerances for high-speed rail track alignment and require periodic realignment maintenance.

4.5. North-South Velocity Field from GNSS

Due to the near-polar orbit geometry of Sentinel-1, InSAR is inherently insensitive to north-south surface displacements. The GNSS network provides direct measurements of the NS velocity component at each station. The GNSS-derived NS velocity field (Figure 15d) reveals important deformation patterns complementing the InSAR-derived EW and vertical fields.

In the subsidence zone, GNSS north-south velocities show variations of several mm/yr correlating with proximity to subsidence centers. For the THSR, which runs in an approximately north-south direction, the NS velocity component represents along-track deformation that could cause differential axial strain along the viaduct. Continuous GNSS stations show NS velocities ranging from −4.2 to +6.1 mm/yr relative to GS26, while campaign GNSS stations show −2.6 to +3.7 mm/yr.

Table 3. Summary of GNSS station categories and velocity statistics.

Type	Count	E range (mm/yr)	N range (mm/yr)	U range (mm/yr)	Reference
Continuous GNSS	44	−1.6 to +5.5	−4.2 to +6.1	−40 to +13	GS26
Campaign GNSS	42	−2.9 to +5.8	−2.6 to +3.7	−37 to +3	GS26
Leveling	922	—	—	−72.5 to +29.0	Class II

Table 3. Summary of GNSS station categories and velocity statistics.

Type	Count	E range (mm/yr)	N range (mm/yr)	U range (mm/yr)	Reference
Continuous GNSS	44	−1.6 to +5.5	−4.2 to +6.1	−40 to +13	GS26
Campaign GNSS	42	−2.9 to +5.8	−2.6 to +3.7	−37 to +3	GS26
Leveling	922	—	—	−72.5 to +29.0	Class II

The combination of InSAR-derived EW velocities with GNSS-derived NS velocities enables a more complete characterization of the horizontal deformation field than either technique alone. This integrated approach is essential for accurately assessing the structural risks to linear infrastructure traversing subsidence zones.

5. Discussion

5.1. Mechanism of Subsidence-Induced Horizontal Displacement

The convergent horizontal displacement pattern observed around the subsidence centers is a direct mechanical consequence of aquifer compaction, a phenomenon that has been progressively recognized and documented in the geomechanics literature over the past three decades. The fundamental theoretical basis was established by Helm [19], who derived analytical solutions for horizontal aquifer movement in a Theis-Thiem confined system, demonstrating that radial groundwater flow toward a pumping well induces horizontal displacement of the aquifer skeleton through drag forces. This work provided the first rigorous framework showing that horizontal displacement is an inherent component of aquifer-system deformation, not merely an incidental byproduct.

Building upon this theoretical foundation, Burbey [14] demonstrated through field observations in Las Vegas Valley that faults in basin-fill deposits significantly amplify horizontal deformation by creating discontinuities in the stress field, resulting in horizontal-to-vertical displacement ratios (uₕ/w) of 0.05 to 0.20. Burbey [20] further showed that neglecting horizontal strain leads to systematic overestimation of specific storage and compaction in confined aquifer systems, because horizontal strain contributes a significant fraction of the total volumetric strain. In a subsequent landmark field experiment, Burbey et al. [21] directly measured three-dimensional deformation during municipal pumping in the Virgin River Valley, Nevada, recording approximately 8 mm of horizontal displacement within 22 days of pumping. That study revealed a slowly migrating outward wave of compressional strain transmitted through an 83-m-thick brittle unsaturated zone, confirming that horizontal deformation at the surface faithfully reflects aquifer compaction at depth. The companion numerical analysis [22] further established that near the pumping well, horizontal strain components (radial and hoop) can be comparable in magnitude to vertical strain, and that ignoring horizontal strain in one-dimensional models systematically overestimates compaction or underestimates storage.

In the physical mechanism, as compressible clay layers (aquitards) undergo vertical compression due to increasing effective stress from groundwater withdrawal, the Poisson effect generates coupled horizontal contraction within the compacting layer, which is transmitted to the surface [14]. In an alluvial aquifer system, the vertical strain (ε_z) induced by effective stress increase is related to the horizontal strain (ε_h) through the Poisson ratio (ν) of the aquitard material. While in a fully laterally unconfined medium the horizontal contraction would be uniform, in practice the subsurface horizontal boundary conditions imposed by surrounding strata create strain gradients that cause the ground surface to converge toward the subsidence center. The horizontal displacement (uₕ) is therefore proportional to the spatial gradient of vertical displacement (∂w/∂r), which explains the key observation in this study and prior work [14,21] that maximum horizontal displacement occurs where the vertical gradient is steepest.

Independent observational evidence from several global case studies corroborates this mechanical coupling. Bawden et al. [16] used InSAR and GPS in the Los Angeles Basin to detect seasonal horizontal oscillations of up to 7 mm associated with groundwater withdrawal and reinjection in the Santa Ana Basin, accompanied by 55 mm of vertical seasonal oscillation and 12 mm/yr of long-term subsidence, yielding an effective uₕ/w ratio of approximately 0.10. Pacheco-Martínez et al. [15] documented significant land subsidence and associated horizontal displacements in the Aguascalientes Valley, Mexico, where excessive groundwater extraction produced an estimated uₕ/w ratio of approximately 0.15 and generated ground fissures attributed to extensional horizontal strain at the margins of the subsidence bowl. In the broader review by Galloway and Burbey [23], horizontal deformation was identified as one of the critical yet under-studied components of aquifer-system compaction, with the authors emphasizing the need for three-dimensional deformation monitoring to fully characterize subsidence hazards. Hernandez-Marin and Burbey [24] used three-dimensional finite-element models to demonstrate that Quaternary faults in Las Vegas Valley create zones of concentrated stress that amplify horizontal displacement and trigger earth fissures where tensional stress from pumping-induced horizontal strain exceeds the tensile strength of near-surface sediments.

The Poisson ratio of the compacting materials is the key parameter governing the horizontal-to-vertical displacement ratio. Laboratory and field studies indicate that saturated clay exhibits Poisson ratios of 0.4–0.5, approaching the incompressible limit, while consolidating or partially drained soft clay ranges from 0.3–0.4, and sandy deposits exhibit lower values of 0.2–0.3 [25,26]. The observed horizontal-to-vertical velocity ratio of 0.05 to 0.15 in the Changhua-Yunlin area is consistent with theoretical predictions for aquitard-dominated compaction with effective Poisson ratios of 0.3–0.4, which aligns well with the known lithological composition of the Choushui River alluvial fan that contains deep compressible clay layers interbedded with sand and gravel aquifers [3,4]. This ratio is comparable to values reported in Las Vegas (0.05–0.20) [14], Los Angeles (approximately 0.10) [16], and Aguascalientes (approximately 0.15) [15], despite differences in geological settings and pumping regimes, suggesting that the Poisson-effect mechanism is a robust and universal feature of pumping-induced aquifer compaction in unconsolidated alluvial systems. Notably, the relatively lower end of the observed ratio in this study area (0.05–0.15 compared with the upper range of 0.20 in Las Vegas) may reflect the greater lateral confinement imposed by the thick and laterally extensive clay aquitards in the Choushui River alluvial fan, which more effectively constrain horizontal deformation relative to the fault-segmented basin-fill setting of Las Vegas [14,24].

It is important to note the distinction between elastic and plastic deformation regimes. When effective stress remains below the preconsolidation stress, compaction is elastic and largely reversible, and the Poisson-effect-driven horizontal displacement is similarly recoverable. However, when effective stress exceeds the preconsolidation stress, as is likely the case in the heavily pumped Changhua-Yunlin area, compaction becomes inelastic and permanent [20,23]. Under such inelastic conditions, the horizontal deformation may become irreversible and potentially more damaging, as the material undergoes permanent grain rearrangement and fabric collapse that can amplify lateral strain beyond elastic predictions. The spatial pattern of horizontal velocities directed radially toward the subsidence centers, as observed in this study, confirms the compaction-driven origin of the horizontal displacement field, as opposed to tectonic deformation which would exhibit a more uniform regional pattern consistent with the known plate convergence direction in Taiwan [27,28].

5.2. Implications for THSR Structural Safety

The implications for THSR structural safety are significant. The rail alignment traverses two major subsidence centers where both vertical and horizontal differential displacements are substantial. For the THSR operating at speeds up to 300 km/h, even millimeter-level track irregularities can affect ride quality and safety.

In the Xizhou zone (Changhua), vertical subsidence rates of 15–25 mm/yr combined with EW velocity gradients of 3–5 mm/yr over 5 km translate, over a 50-year design life, to cumulative differential vertical settlement of 0.75–1.25 m and horizontal displacement of 15–25 cm. In the Tuku zone (Yunlin), vertical subsidence exceeds 50 mm/yr with EW velocity differences of up to 8 mm/yr, implying cumulative horizontal displacements of up to 40 cm over 50 years.

These findings underscore the need for: (1) continuous geodetic monitoring using integrated InSAR, GNSS, and leveling; (2) structural assessment considering both vertical and horizontal components; and (3) adaptive maintenance strategies accounting for time-dependent accumulation of three-dimensional strain along the THSR corridor.

5.3. Comparison of EW and Vertical Calibration Strategies

An important finding of this study is the fundamentally different calibration strategies required for the EW and vertical components, and the contrasting degrees of improvement achieved by each. Understanding why the EW correction yields a dramatically higher final accuracy (R² = 0.992, RMSE = 0.24 mm/yr) than the vertical correction (R² = 0.898, RMSE = 4.85 mm/yr) despite the vertical component having far more control points (922 leveling benchmarks versus 89 GNSS stations) requires consideration of three interacting factors: the nature of the systematic errors, the correction methodology, and the inherent accuracy limits of each observing system.

First, the contrasting initial accuracies reflect the fundamentally different sensitivity of InSAR to EW and vertical motion. In the 2.5D decomposition geometry of Sentinel-1 (incidence angles of ~34–40°), the vertical component dominates the LOS signal (sensitivity coefficient cos(θ) ≈ 0.77–0.87), whereas the EW component contributes far less (sin(θ) ≈ 0.49–0.64). Consequently, the raw InSAR vertical field already captures the dominant deformation signal (R² = 0.889 against leveling), whereas the EW field is almost entirely obscured by systematic biases (R² = 0.147), because even small long-wavelength errors—on the order of several mm/yr from orbital inaccuracies, ionospheric gradients, and tropospheric ramps [29]—overwhelm the weaker EW signal component.

Second, the correction methodology itself explains much of the performance gap. The EW field is corrected by Ordinary Kriging interpolation of GNSS residuals, a method that by construction forces the corrected field to match the GNSS control points exactly. The resulting near-perfect agreement (R² = 0.992) therefore partly reflects statistical self-consistency: the same GNSS stations serve as both calibration inputs and the primary validation reference. In contrast, the vertical field is corrected by a 2nd-degree polynomial trend surface fitted to the combined leveling and GNSS observations. This approach is deliberately conservative; it removes only the smooth, long-wavelength bias without imposing local constraints at individual benchmarks and is validated against an independent leveling dataset not used in deriving the polynomial coefficients. The more modest improvement (RMSE reduction from 6.03 to 4.85 mm/yr, 17%) therefore provides a more honest estimate of absolute accuracy. Furthermore, the effectiveness of the Kriging correction for the EW component depends critically on the spatial density of GNSS stations: as demonstrated by the comparison between CGPS-only (R² = 0.815) and combined CGPS plus campaign GNSS (R² = 0.992) corrections, denser station coverage resolves shorter-wavelength error components that a sparse network cannot constrain [30].

Third, the residual RMSE of 4.85 mm/yr in the corrected vertical field approaches an irreducible noise floor that cannot be eliminated by any external calibration strategy. This floor arises from at least three independent sources: (1) atmospheric phase noise from tropospheric turbulence, which contributes 1–3 mm/yr to InSAR time-series velocity estimates and cannot be fully removed without high-resolution meteorological corrections [31]; (2) spatial scale mismatch between the 20-m InSAR pixel and point-support leveling benchmarks, which introduces representativeness errors of 1–3 mm/yr in areas with rapid spatial gradients in subsidence; and (3) temporal sampling differences between the annual leveling snapshots and the continuous multi-year InSAR velocity estimation, which introduces additional uncertainty of 1–4 mm/yr when subsidence rates vary seasonally or inter-annually. Crucially, the inherently lower vertical accuracy of GNSS itself—with vertical errors typically 1.5–2 times larger than horizontal errors due to the asymmetric satellite geometry above the observer’s horizon (VDOP/HDOP ≈ 1.5)—limits the precision of the polynomial trend surface fit and prevents GNSS vertical observations from serving as a high-accuracy anchor in the same way that GNSS horizontal velocities anchor the Kriging EW correction. In summary, the superior numerical performance of the EW calibration reflects a combination of methodological design (Kriging vs. polynomial), the self-referencing nature of Kriging validation, and the availability of a dense, high-accuracy GNSS horizontal velocity field. The modest improvement in the vertical field, by contrast, represents a genuine physical limit imposed by atmospheric noise, scale mismatch, and GNSS vertical precision constraints. The leveling-based RMSE of 4.85 mm/yr should therefore be interpreted not as a calibration shortcoming, but as a realistic and externally validated accuracy estimate for InSAR vertical velocity fields in low-relief alluvial fan settings with complex spatiotemporal subsidence patterns.

5.4. Limitations and Future Work

A fundamental limitation is the insensitivity of InSAR to north-south displacements, which means that the horizontal velocity field derived from InSAR alone captures only the east-west component. Although GNSS provides NS velocities at discrete stations, the spatial resolution of the NS field is limited by the GNSS network density. Future studies incorporating ascending and descending data from multiple satellite platforms (e.g., ALOS-2, COSMO-SkyMed) with different orbit geometries could potentially improve the sensitivity to the NS component.

The campaign GNSS data showed degraded performance after vertical correction (RMSE increased from 5.60 to 7.70 mm/yr), likely due to short observation periods (1–3 days) and incomplete temporal overlap with the InSAR period. Future work should prioritize campaign GNSS observations with longer occupation times and temporal alignment with the InSAR processing window.

Additionally, the temporal coverage represents only a snapshot of the ongoing subsidence process. Long-term monitoring with continuous InSAR processing is essential for capturing temporal variations in subsidence rate related to groundwater management policies and climate variability.

6. Conclusions

This study presents an integrated geodetic framework combining MT-InSAR, GNSS, and precise leveling to characterize the full three-dimensional (3D) surface deformation field in the Choushui River alluvial fan. Beyond providing high-precision monitoring, this research elucidates the mechanical coupling between vertical compaction and lateral strain, offering critical insights for infrastructure resilience. The principal conclusions are as follows:

• Methodological Breakthrough in Horizontal Velocity Calibration: By implementing an Ordinary Kriging calibration using combined continuous and campaign GNSS data, the InSAR-derived EW velocity field accuracy was dramatically improved, with R² increasing from 0.147 to 0.992 and RMSE decreasing from 4.24 to 0.24 mm/yr. This demonstrates that even when horizontal signals are initially obscured by systematic biases, dense ground-truth constraints can recover millimeter-level tectonic and anthropogenic signals.
• Comprehensive 3D Characterization: The calibrated vertical velocity field revealed maximum subsidence exceeding 60 mm/yr in the Tuku area. By integrating GNSS-derived north-south velocities (ranging from −4.2 to +6.1 mm/yr), this study overcomes the inherent InSAR “blind direction,” achieving a complete 3D deformation profile essential for linear infrastructure monitoring.
• Mechanical Evidence of Aquifer Compaction: The identification of a radially convergent horizontal velocity pattern (2–10 mm/yr) pointing toward subsidence centers provides direct observational evidence of the Poisson effect during aquitard compaction. The observed horizontal-to-vertical displacement ratios (0.05–0.15) align with theoretical mechanical models for clay-dominated aquifer systems.
• Quantification of Long-term Engineering Risks: A critical finding is that the maximum lateral hazard to the THSR occurs at the inflection points of the subsidence bowl—where the vertical gradient is steepest—rather than at the subsidence center itself. The horizontal strain rate in these transition zones reaches approximately 10⁻⁶/yr.
• Implications for High-Speed Rail Maintenance: Projections indicate that the THSR will experience 15–40 cm of cumulative horizontal displacement over its 50-year design life. These results underscore the necessity for structural health monitoring strategies to evolve beyond a traditional vertical-subsidence focus and incorporate 3D strain accumulation into maintenance planning and track realignment.

The integrated approach established in this study provides a robust, transferable framework for monitoring land subsidence in alluvial fans worldwide, offering a scientific basis for both groundwater management and the long-term safety assessment of major civil infrastructure.