Spatiotemporal evolution of land subsidence in coastal reclamation areas: coupling anthropogenic land-use intensity and delayed soft-soil compression

1 Introduction

Rapid global urbanization and intensive coastal-land development have made land subsidence a major geohazard for coastal metropolises (Wu et al., 2022; Hadianti and Dewanto, 2023; Huang L. et al., 2025). The risk is particularly acute in estuarine and coastal reclamation areas, where artificial reclamation, complex marine strata, and thick soft-soil deposits create a highly vulnerable deformation environment (Chen, 2011; Jiang et al., 2020; Van Maren et al., 2025). Unlike inland subsidence dominated by aquifer depletion, deformation in newly reclaimed areas is persistent, spatially heterogeneous, and nonlinear. It reflects the primary consolidation of dredged fills, dissipation of excess pore-water pressure, and secondary consolidation of alluvial and marine clays.

Process-based prediction of reclamation-induced settlement has traditionally been grounded in effective stress consolidation theory and geotechnical constitutive modeling. Classical studies established the importance of delayed consolidation and the relationship between time-dependent and stress-dependent compressibility in marine clays (Bjerrum, 1967; Mesri and Godlewski, 1977). Studies of Chinese coastal soft soils further show that secondary consolidation is affected by loading level, stress history, compression characteristics, water content, and plasticity (Chen, 2011; Li et al., 2012; Qin et al., 2020; Liu et al., 2022).

In the present study, secondary consolidation is therefore regarded as part of the physical background of long-term soft-soil deformation rather than as a process directly resolved by the numerical constitutive formulation. The regional model developed below uses a site-conditioned equivalent representation of delayed compression and does not explicitly reproduce material-point creep evolution or microstructural changes.

In parallel, time-series InSAR has become an important tool for detecting regional subsidence and investigating its potential drivers. The permanent-scatterer method and its subsequent developments provide the methodological basis for long-term deformation monitoring using coherent radar targets (Ferretti et al., 2001; Crosetto et al., 2016). Global observations have demonstrated the widespread occurrence of rapid subsidence in low-lying coastal cities (Wu et al., 2022). Regional studies have further related subsidence patterns to Quaternary sediment thickness, reclaimed land, built-up areas, groundwater conditions, and engineering activities (Ren et al., 2024; Zhou et al., 2025). Persistent scatterer InSAR (PS-InSAR) observations have also been combined with geological, hydrological, land-use, and infrastructure-related variables in statistical and machine-learning frameworks for subsidence susceptibility mapping and regional risk assessment (Qiao et al., 2024; Zhao et al., 2024). Recent studies have further integrated Sentinel-1 observations with ground-based measurements and data-driven models to support regional subsidence assessment (Zhao et al., 2024; Ma et al., 2025).

These studies have substantially improved the monitoring, attribution, susceptibility assessment, and statistical prediction of land subsidence. However, most InSAR-based multi-factor approaches treat geology, land use, groundwater, and infrastructure as explanatory or predictive variables without explicitly translating the evolution of anthropogenic development into geomechanical loading conditions (Zhang et al., 2023; Qiao et al., 2024; Zhao et al., 2024). Conversely, process-based geotechnical models can simulate stress transfer and soil compression but often simplify construction loading as static or spatially averaged boundary conditions. The mechanical interaction between thick underconsolidated marine soft soils and progressively intensified land use therefore remains insufficiently represented in regional-scale prediction.

Land-Use Intensity (LUI) provides a potential link between urban spatial evolution and regional geomechanical loading. Changes from sparsely developed reclaimed land to logistics areas, commercial districts, storage yards, and port infrastructure are accompanied by increasing fill, structural, operational, and construction-related loads. However, land-use indicators are commonly employed as categorical correlates or susceptibility factors rather than being converted into spatially distributed and temporally evolving surface stresses (Hadianti and Dewanto, 2023; Qiao et al., 2024; Zhao et al., 2024). Establishing such a conversion is particularly important in newly reclaimed areas, where incremental anthropogenic loading is superimposed on marine soft soils that may still be undergoing self-weight consolidation and delayed compression.

The Ningbo coastal zone provides a representative setting for examining this coupled process. Recent InSAR investigations indicate that, following restrictions on groundwater extraction, active subsidence has increasingly concentrated in coastal reclamation and intensive development zones (Qiao et al., 2024; Zhou et al., 2025; Ma et al., 2025). The Meishan Bonded Port Area is characterized by spatially variable, thick underconsolidated marine deposits and rapidly evolving port and urban land use. In 2020, localized annual subsidence rates in the southern reclamation sector exceeded 40 mm, posing potential risks to terminals, storage facilities, underground infrastructure, and drainage systems. These geological and development conditions require a framework capable of representing both subsurface compressibility and the spatial evolution of anthropogenic loading.

Accordingly, this study develops an integrated regional prediction framework that combines a borehole-constrained three-dimensional geological and finite-element model, multi-level settlement monitoring, PS-InSAR observations, and LUI-based surface-load inversion. The contribution of the framework lies in three aspects. First, the spatial variation of the principal soft- and stiff-soil units is retained in the full-domain three-dimensional model, allowing geological susceptibility to be represented more explicitly than in uniform-layer or two-dimensional regional approximations. Second, a layer-specific equivalent parameterization is used to represent the first-order effect of delayed compression in underconsolidated marine soft soils under regional modeling constraints. Third, categorical LUI is translated into spatially distributed and temporally evolving equivalent surface loads, thereby introducing anthropogenic development directly into the geomechanical boundary conditions. The framework is evaluated against stratified settlement observations and regional PS-InSAR patterns and is subsequently applied to investigate the spatiotemporal evolution and combined geological–anthropogenic controls of subsidence from 2025 to 2035. The overall research workflow is summarized in Figure 1.

2 Materials and methods

2.1 Geomechanical discretization and multi-level stratigraphic evaluation

2.1.1 Principles of stratigraphic generalization

The construction of a regional 3D geological model provided the spatial framework for the subsequent geomechanical analysis. In a heterogeneous reclamation environment, stratigraphic generalization must preserve the principal geological boundaries and hydro-mechanical contrasts while maintaining a level of complexity suitable for regional-scale numerical calculation. Multi-source borehole data were integrated to classify the strata according to geological age, genesis, soil type, and consolidation state. Water content, permeability, and mechanical behavior were also considered (Zhong et al., 2006; Chen et al., 2017; Zhang et al., 2023; Ji D. et al., 2024; Pusacker et al., 2024).

Local sub-layers with comparable geological and mechanical characteristics were grouped into generalized stratigraphic units. The top and bottom elevations of these units were determined from the selected boreholes and spatially interpolated to construct continuous three-dimensional stratigraphic surfaces. Consequently, the resulting model preserved the lateral variation in layer elevation and thickness across the study area rather than applying a single uniform vertical profile to the entire domain.

The generalization procedure emphasized the continuity of borehole correlations, consistency among geological cross-sections, and preservation of the principal contrasts between soft- and stiff-soil units (Ji G. et al., 2024; Zhang et al., 2024). In particular, the spatial variation of the thick underconsolidated soft-soil layers was retained because it constitutes one of the principal geological controls on regional settlement. This treatment provided a balance between local geological detail and the computational requirements of the full-domain three-dimensional finite-element model.

2.1.2 Evaluation using high-resolution multi-level settlement monitoring

To evaluate the plausibility of the geomechanical discretization and the assigned baseline constitutive parameters before regional upscaling, a localized evaluation was conducted using a high-resolution multi-level settlement monitoring array in the Meishan Bonded Port Area. Unlike conventional surface-leveling observations, which integrate deformation over the entire soil column, the multi-level monitoring system provides depth-specific compression records and therefore enables the deformation behavior of individual stratigraphic intervals to be examined.

A spatially co-located borehole, HMG06 (Figure 2), was selected to define the stratigraphic sequence of the local vertical model constructed in PLAXIS. The associated multi-level settlement records covered 2018-2021 and comprised three monitored depth intervals: 1.5–5.3 m, 5.3–17.5 m, and below 17.5 m. Historical overburden loading and seasonal groundwater fluctuations were imposed, and the simulated annual compression of the three intervals was compared with the corresponding measurements. The HMG06 observations were used to evaluate the combined stratigraphic and constitutive representation after the equivalent correction scheme had been specified; they were not used in a sensitivity sweep or inverse optimization of the correction multipliers. The HMG06 record spans 2018-2021, whereas the regional PS-InSAR product used in Section 3.2 represents 2022.

Descriptive evaluation used annual total-compression relative errors and the 12 repeated layer-year pairs to calculate mean absolute error (MAE), root mean square error (RMSE), and mean bias. Because the observations came from one location over repeated years, they were not treated as independent replicates; inferential hypothesis tests and confidence intervals based on independent samples were therefore not used.

2.1.3 PS-InSAR data and processing

The regional deformation information used for comparison was the processed subsidence-rate product for 2022. It was derived from a RADARSAT-2 C-band Extra Fine VV archive acquired in descending, right-looking geometry. The archive comprised approximately 90 scenes acquired between 9 December 2016 and 16 September 2023; the basic repeat interval was about 24 days and increased to approximately 48 days where acquisitions were absent. Incidence angles ranged from approximately 21°–33°, with 27.8° used as the representative angle for LOS-to-vertical conversion. GAMMA was used for precise co-registration, DEM-related phase removal, and differential-interferogram generation. StaMPS was then used for persistent-scatterer screening and time-series inversion under a conventional small-baseline framework. Temporal and perpendicular baselines were limited to 72 days and 200 m, respectively.

Persistent-scatterer candidates were identified primarily from backscatter intensity and coherence, with amplitude-dispersion and phase-stability information used in subsequent screening. A 32-pixel Goldstein filter and SNAPHU phase unwrapping were applied, and the retained point density exceeded 500 PS/km². No external atmospheric-correction product was applied; residual long-wavelength phase components were addressed within the standard time-series filtering procedure. A phase-stable reference area with dense coherent targets and separation from the principal subsidence center was used to define relative deformation. The processing yielded LOS displacement and deformation-rate results. Assuming horizontal motion was small relative to vertical deformation, the approximate vertical rate was calculated as VU = VLOS/cos θ. Negative vertical rates denote subsidence. PS-InSAR velocity precision was characterized using the per-point standard deviation of the mean LOS deformation rate supplied by StaMPS. Across the study-area PS points, the median LOS-rate standard deviation was approximately 2.6 mm/year. After LOS-to-vertical projection using the representative incidence angle, the corresponding median vertical-rate standard deviation was approximately 3.0 mm/year. The resulting 2022 approximate vertical subsidence-rate product was used for the regional comparison. Table 1 summarizes the PS-InSAR dataset and key processing parameters.

TABLE 1. Summary of PS-InSAR data and key processing parameters.

PS-InSAR, velocity precision was evaluated using the StaMPS-derived per-point standard deviation of the mean deformation rate.

2.2 Multi-factor predictive framework formulation

2.2.1 Constitutive modeling and site-conditioned equivalent parameterization of delayed soft-soil compression

Regional settlement was simulated using a full-domain three-dimensional finite-element model established in PLAXIS 3D. The spatial geometry of the model was derived from the generalized stratigraphic framework described in Section 2.1, thereby preserving the lateral variation in the thickness and distribution of the principal soft- and stiff-soil units. The model was discretized using the automatic mesh-generation function in PLAXIS 3D, producing an unstructured tetrahedral mesh conforming to the three-dimensional stratigraphic interfaces.

The regional geometry was controlled by 17 selected boreholes. Boreholes were retained at an approximate 1,000 m spacing where coverage permitted, and the top and bottom elevations of each generalized unit were interpolated between the borehole controls. The horizontal model extent covered approximately 36.5 km² and followed the irregular Meishan study-area boundary shown in Figures 3–6, while the model base followed the lower boundary of the fourth stiff-soil unit; consequently, the computational domain was a stratigraphically varying three-dimensional volume rather than a constant-depth rectangular block.

The ground surface was partitioned into 500 m × 500 m load cells, and the LUI-derived equivalent loads were applied vertically in the negative Z direction. Automatic tetrahedral meshing was followed by local refinement at load-cell transitions and where stratigraphic surfaces changed abruptly. The base was fixed in X, Y, and Z; lateral boundaries were restrained normal to the boundary in the horizontal directions and remained free vertically. Initial stresses were generated with the K0 procedure, with K0,eq = 0.80K0,base only in the designated shallow underconsolidated units. Consolidation phases used the assigned unit permeabilities and a time-dependent phreatic boundary; the Meishan groundwater depth varied sinusoidally from 1.25 to 1.75 m below ground surface over a 1-year cycle. The staged calculation comprised the 2022 baseline state and forward phases ending in 2025, 2030, and 2035, with the corresponding LUI-load field updated at each target year. Mesh density was generated automatically from the irregular stratigraphic geometry and the specified local-refinement rules.

Long-term settlement in coastal reclamation areas includes primary consolidation during excess pore-water-pressure dissipation and secondary consolidation or creep under sustained effective stress (Bjerrum, 1967; Mesri and Godlewski, 1977; Yin and Graham, 1989). Laboratory and field studies further show that secondary-compression behavior varies with loading level, stress history, water content, compressibility, and soil structure (Qin et al., 2020; Liu et al., 2022). Recent work also demonstrates that soft-soil deformation may involve coupled creep and nonlinear consolidation with evolving material properties (Ouria, 2024). Explicit creep formulations, such as the Soft Soil Creep model, require layer-specific long-term parameters (Vermeer and Neher, 1999). Given the regional scale and the need for a consistently constrained formulation across the full three-dimensional domain, the conventional Mohr-Coulomb model was retained and a fixed, layer-level equivalent parameter scheme was used to approximate the first-order magnitude of delayed compression. The M-C law itself was not converted into a creep constitutive law, and the scheme does not reproduce material-point creep evolution.

The estimated consolidation degree and the parameter-reduction factors represent distinct quantities. The 70%–80% value is a site-conditioned engineering screening range adopted from the depositional setting, borehole-based consolidation classification, and persistence of shallow compression; it is not presented as a direct site-wide measurement. Experimental and field studies show that underconsolidated marine clay retains excess pore pressure and that its undrained strength evolves with degree of consolidation (Towhata and Kim, 1990; Karakouzian et al., 2003). The range therefore provides a physically plausible descriptor of incomplete consolidation, but not a universal threshold. For only the identified underconsolidated units, the equivalent values were defined as K_0,eq = 0.80K_0,base, E_eq = (0.70–0.80)E_base, and φ_eq = (0.70–0.80)φ_base. Thus, the operational elastic modulus and internal friction angle were reduced by 20%–30%, not by 70%–80%.

After these reductions, apparent cohesion was set to c _eq = (1.20–1.30)c _base to suppress excessive localized shear yielding and maintain numerical stability in the regional calculation. Accordingly, it is treated solely as a numerical stabilizer within the surrogate scheme rather than as a measured increase in intrinsic soil cohesion.

The model inputs draw on distinct lines of evidence. The multi-source borehole and site-investigation dataset supports the identification of the targeted shallow underconsolidated units and the baseline properties reported in Table 2; it does not independently prescribe the approximate consolidation-state descriptor or the exact M-C multipliers. The cited literature establishes the physical importance of stress- and time-dependent compression, incomplete consolidation, and evolving soft-soil behavior (Bjerrum, 1967; Mesri and Godlewski, 1977; Yin and Graham, 1989; Towhata and Kim, 1990; Karakouzian et al., 2003; Qin et al., 2020; Liu et al., 2022; Ouria, 2024), but none of these sources prescribes the present factors as a universal relationship. The exact factors are therefore treated as site-conditioned engineering assumptions for regional screening, selected with reference to the local consolidation classification, baseline parameter ranges, and stable numerical behavior. Their role is to define an equivalent regional approximation rather than measured constitutive parameters.

TABLE 2. Borehole-derived parameter envelopes, recommended shear-strength values, and operational inputs for the simplified strata in Meishan Bonded Port Area.

w, γ, e, G_s, and E_s are borehole-derived spatial envelopes; C_k and φ_k are recommended standard values. Baseline inputs vary by borehole-controlled generalized unit. In designated shallow underconsolidated first-soft-soil cells, the declared correction factors give Es = 1.20–3.44 MPa, c = 19.8–21.5 kPa, φ = 8.7°–9.9°, and K_0,eq = 0.80K_0,base; other cells retain baseline inputs.

The same layer-level correction rules were applied to all identified shallow underconsolidated units in the regional model without location-specific retuning, whereas baseline values were retained for stiff or non-target strata. Evaluation therefore has two complementary levels: the HMG06 record tests depth-resolved compression at one monitored column, and the 2022 regional PS-InSAR comparison tests the spatial pattern over 3,872 grid cells. The latter yielded a quadratic weighted Cohen kappa of 0.684, an intersection-over-union of 0.510 for the high-subsidence zone, and a centroid offset of approximately 0.81 km. Together, these comparisons provide complementary depth-resolved and regional spatial assessments, while model-form uncertainty associated with the M-C surrogate remains.

Transferability therefore resides in the workflow - borehole-constrained three-dimensional geometry, explicit declaration of equivalent factors, depth-resolved local evaluation, and a regional spatial consistency assessment - rather than in the numerical multipliers themselves. The present factors are restricted to regional screening of recently reclaimed areas with underconsolidated marine soft soils and geological, groundwater, loading, and consolidation conditions comparable to those represented in Meishan. Applications to other sites require recalibration against local laboratory or field observations and independent validation; materially different soft-soil thickness, consolidation state, groundwater regime, or loading intensity precludes direct use of the present factors.

2.2.2 Borehole selection and stratigraphic parameterization

The regional model used 17 selected boreholes, although their spatial density and survey dates vary. Boreholes drilled after 2018 were prioritized in densely investigated areas. Where multiple recent boreholes were available, representative boreholes were selected at intervals of approximately 1,000 m. In areas without post-2018 boreholes, the most recent available boreholes were retained to maintain spatial coverage. Figure 3 shows the final borehole distribution used in the analysis.

According to the latest stratigraphic classification guidelines of Ningbo City, the strata were grouped into the First Soft Soil Layer, First Stiff Soil Layer, Second Stiff Soil Layer, Second Soft Soil Layer, Third Stiff Soil Layer, Third Soft Soil Layer, and Fourth Stiff Soil Layer, as summarized in Table 3. The original borehole layers were reassigned to these generalized units according to their stratigraphic numbers and geological characteristics.

TABLE 3. Names and layer numbers of stratigraphic divisions.

The elevations of the upper and lower boundaries of each generalized unit were retained at the borehole locations and subsequently interpolated to generate spatially continuous stratigraphic surfaces. Therefore, the burial depth and thickness of each unit varied across the three-dimensional model according to the borehole-controlled geological geometry. The physical and mechanical properties of the original sub-layers assigned to the same generalized unit were statistically summarized, and representative parameter ranges or recommended values were adopted for the numerical analysis. These values are presented in Table 2 and served as the material-property inputs for the regional finite-element model. Table 2 distinguishes the roles of the reported values: water content, unit weight, void ratio, specific gravity, and compression modulus are borehole-derived spatial envelopes; cohesion and friction angle are recommended standard values. Baseline inputs were assigned locally by borehole-controlled generalized unit rather than as one domain-wide value. In cells designated as shallow underconsolidated first soft soil, the operational envelopes resulting from the declared factors were E_s = 1.20–3.44 MPa, c = 19.8–21.5 kPa, and φ = 8.7°–9.9°, with K_0,eq = 0.80K_0,base; non-target cells retained the baseline inputs.

2.2.3 Geospatial inversion of anthropogenic surface loads based on land-use intensity (LUI)

Based on the official Main Data Bulletin of the Third National Land Survey of Ningbo City (Ningbo Leading Group Office for the Third National Land SurveyNingbo Municipal Bureau of Natural Resources and Planningand Ningbo Municipal Bureau of Statistics, 2021), construction land proportion within the Meishan Bonded Port Area was classified into four LUI grades (Table 4; Figure 4). Grade 1 (<20% construction land) represents sparsely developed or buffer zones, whereas Grade 4 (>60% construction land) denotes areas of intensive urbanization, including commercial clusters, logistics hubs, and port infrastructure.

TABLE 4. Correspondence between construction land density and grade.

To convert categorical LUI into a geomechanical boundary condition, the construction-land proportion grade was first extracted at the borehole locations used in the regional model. For the Meishan area, one representative borehole profile was selected within each of the four LUI grades to provide a grade-specific stratigraphic column rather than a spatially independent validation point. Selection emphasized consistency with the corresponding land-use grade, representativeness of the local generalized stratigraphy, and usable InSAR coverage. For each grade-specific profile, the 2022 PS-InSAR-derived vertical subsidence-rate magnitude was used as the calibration target. The target rates for LUI Grades 1-4 were 4.51, 13.44, 22.63, and 31.97 mm/year, respectively. Starting from a trial surface load, the load was increased when the modeled subsidence-rate magnitude was below the corresponding target and decreased when it was above the target. The retained load was the trial value that minimized the absolute difference between the modeled and PS-InSAR-derived annual rates. The four resulting Meishan calibration points were (G, L) = (1, 1.4), (2, 2.2), (3, 3.8), and (4, 4.8), where G is the LUI grade and L is the back-calculated equivalent load in kPa. An origin-constrained least-squares fit gave L = 1.213G with R ² = 0.984, calculated against the centered total sum of squares; rounding the slope to one decimal place for regional application yielded the Meishan-specific relationship in Equation 1, L = 1.2G. The zero intercept is a boundary constraint indicating that a hypothetical state without an LUI-related development increment carries no corresponding incremental load; it lies outside the observed Grade 1-4 range. The reported R ² describes the linear consistency of the four grade-specific calibration points and is not interpreted as independent predictive validation or as a coordinate-based spatial regression.where L represents the equivalent overlying load in kPa and G is the construction land density grade, ranging from 1 to 4.

L = 1.2 × G (1)

This fitted relationship provides a Meishan-specific approximate conversion from LUI to the surface-load boundary used in the regional model. The resulting 1.2–4.8 kPa range represents settlement-equivalent, area-averaged incremental loads applied to 500 m × 500 m cells. These values are not foundation pressures or the structural load of an individual building. Instead, they aggregate the net regional influence of shallow excavation disturbance, fill, structures, storage and port operations, and construction activity after spatial averaging and settlement-based back-calculation. Their interpretation and use are therefore restricted to the Meishan regional screening model.

Future LUI evolution was represented using a baseline planning scenario. According to the Master Plan for the Meishan Bonded Port City of Ningbo (2010–2030), the planned construction land area will reach approximately 19.99 km² by 2030, corresponding to 52.2% of the total Meishan Island area. The construction land proportion was 32.2% in 2019 according to the Third National Land Survey bulletin. Linear interpolation between the 2019 condition and the 2030 planning target yields an average annual LUI growth rate of 1.82%.

For the period from 2030 to 2035, the same annual rate was extended solely as a baseline scenario assumption to maintain a continuous anthropogenic-loading trajectory. This extension should not be interpreted as an approved post-2030 planning target or as a deterministic forecast of actual development. Slower or accelerated development relative to this trajectory would correspondingly decrease or increase the imposed anthropogenic loading and the predicted local settlement.

2.2.4 Hydrogeological boundary conditions

Groundwater-level variation affects pore-water pressure and the effective stress of shallow saturated soils and therefore influences the ongoing consolidation of the reclaimed deposits (Liu et al., 2019; Gao, 2021). Based on the available geotechnical investigation reports and regional pore-water monitoring information, the baseline groundwater table was set at 1.5 m below ground surface, with a seasonal fluctuation amplitude of 0.25 m. An annual sinusoidal boundary varying between 1.25 and 1.75 m below ground surface was applied in the model.

This boundary represents the dominant seasonal fluctuation observed under the available regional monitoring conditions. Long-term changes associated with future excavation, drainage, reclamation, or groundwater-management activities were not prescribed separately and therefore remain a source of uncertainty in the long-term prediction.

3 Results

3.1 Depth-specific model evaluation at HMG06

Tables 5–7 and Figure 7 present the measured and simulated compression values. The measured and simulated annual total compression values were 13.20 and 13.49 mm in 2018, 10.30 and 9.55 mm in 2019, 17.30 and 15.95 mm in 2020, and 11.80 and 10.89 mm in 2021, respectively. The corresponding relative errors were 2.20%, 7.28%, 7.80%, and 7.71%.

TABLE 5. Descriptive comparison of measured and simulated soil-layer compression. Panel A. Measured and simulated compression at different soil-layer depths from 2018 to 2021.

TABLE 6. Panel B. Annual total compression and relative errors.

TABLE 7. Panel C. Descriptive error metrics based on the 12 repeated layer-year observations.

Across the 12 repeated layer-year observations, MAE was 1.22 mm, RMSE was 1.62 mm, and mean bias, defined as the simulated value minus the measured value, was −0.23 mm. These descriptive metrics indicate agreement at HMG06 but do not constitute independent validation of the full regional domain.

3.2 Regional spatial consistency assessment using PS-InSAR observations

The simulated subsidence-rate distribution for 2022 was compared with the processed regional PS-InSAR-derived subsidence-rate product for 2022 (Figure 5). Both datasets show a principal high-subsidence area in the southern reclamation zone. The comparison was conducted at the delivered spatial support of the processed product. Because the same processed 2022 PS-InSAR product also informed the four-point LUI-load calibration, the 2022 comparison is conservatively treated as a regional spatial-consistency assessment rather than independent validation. The HMG06 depth-specific evaluation uses the separate 2018–2021 monitoring series described in Section 2.1.2.

To quantitatively assess the spatial consistency between the simulated and PS-InSAR-derived results, both maps were converted to a common 100 m evaluation grid. Each grid cell was assigned the subsidence-rate class occupying the largest proportion of its area, and 3,872 valid cells within the common study area were retained for comparison. As summarized in Table 8, the exact class agreement was 23.32%, while 68.70% and 90.19% of the evaluation cells differed by no more than one and two legend classes, respectively. The signed cell-by-cell class difference was calculated on the common grid, with a mean absolute value of 1.21 classes, and the quadratic weighted Cohen’s kappa coefficient was 0.684. For the high-subsidence zone defined as ≥20 mm/year, the intersection-over-union was 0.510, and the centroid offset between the simulated and PS-InSAR-derived zones was approximately 0.81 km. These results indicate that most discrepancies were concentrated in adjacent subsidence classes and that the model reproduced the principal southern high-subsidence center and the overall regional gradient, although local class-level differences remained. The 2022 comparison evaluates the baseline spatial pattern of the model, whereas the 2025–2035 results represent scenario-based projections under the prescribed geological, groundwater, and anthropogenic-loading conditions. Because the two datasets have different observational supports, the class-based common-grid comparison characterizes regional consistency without implying pointwise millimeter-scale equivalence.

TABLE 8. Spatial classification consistency between the simulated and PS-InSAR-derived subsidence-rate maps in 2022.

Both maps were converted to a common 100 m grid, and each grid cell was assigned the subsidence-rate class occupying the largest proportion of its area. The high-subsidence zone was defined using a threshold of 20 mm/year.

3.3 Long-term spatiotemporal predictions and temporal attenuation (2025–2035)

Forward simulations were performed to project annual subsidence rates and cumulative deformation for 2025, 2030, and 2035. The resulting spatial distributions are shown in Figure 6, and the areal proportions of the annual and cumulative subsidence classes are summarized in Tables 9, 10 and Figures 8, 9, respectively.

TABLE 9. Area proportion changes of different subsidence rates (2025–2035).

TABLE 10. Area proportion variation trends of cumulative subsidence (2025–2035).

3.4 Analysis of subsidence variation trends

A defining characteristic of the baseline-scenario projection is the persistent spatial entrenchment of the primary subsidence funnel in the southern reclamation zone, which is linked to the distribution of exceptionally thick underconsolidated marine clays. Under the prescribed assumptions, modeled subsidence intensity attenuates over time. Table 9 indicates that the maximum annual rate at the core of the funnel is projected to be 35–40 mm/year in 2025 and to decrease to 15–20 mm/year by 2035. At the same time, high-velocity zones contract markedly: areas exceeding 15 mm/year shrink, and zones exceeding 25 mm/year disappear by 2035, while low residual subsidence zones of 5–10 mm/year expand to more than 41.3% of the study area. This modeled deceleration is compatible with delayed compression but may also reflect the prescribed staging, loading, groundwater, drainage, and boundary conditions; the equivalent M-C formulation does not constitutively resolve creep evolution. Additional loading associated with LUI expansion may induce localized primary consolidation, while ongoing excess-pore-water-pressure dissipation and delayed compression may contribute to the longer-term deformation pattern. Despite the reduction in annual rates, cumulative strain remains severe. By 2035, the peak cumulative subsidence at the southern center reaches 400–450 mm, and the difference between the southern center and the relatively stable northern sector exceeds 400 mm.

The results of cumulative land subsidence, as shown in Table 10 and Figure 9, indicate that the southern subsidence center exhibits a typical “subsidence funnel” pattern. In 2025, the maximum cumulative subsidence is between 200 and 250 mm, with 85% of the area experiencing less than 150 mm. By 2035, the maximum cumulative subsidence increases to 400–450 mm, although 80% of the area remains below 300 mm. Notably, regional subsidence disparity is significant, with the difference between the maximum and minimum cumulative subsidence exceeding 400 mm in 2035. In terms of spatial distribution, the proportion of the area with cumulative subsidence below 150 mm decreases from 85.3% in 2025 to 32.7% in 2035, while the proportion with 150–300 mm rises from 14.7% to 46.4%, and areas with over 300 mm account for 20.9%. These results indicate that over time, the extent of areas experiencing greater cumulative subsidence gradually expands, although the overall spatial pattern of higher subsidence in the south and lower in the north remains unchanged.

3.5 Geological and anthropogenic controls on subsidence

Regional subsidence reflects the combined influence of geological susceptibility and anthropogenic loading. The predicted subsidence pattern, regional LUI distribution, and stratified compression at the representative southern subsidence center (Figures 4, 10, 11) were jointly examined to interpret the respective roles of these controls. The analysis is interpretive and does not quantitatively separate geological and anthropogenic contributions or identify a causal interaction term.

3.5.1 Natural geological control: compression of underconsolidated soft soils

From a geological perspective, the distribution, thickness, and initial consolidation state of deep marine soft soils define the baseline susceptibility to land subsidence. At the representative southern subsidence center, the stratum-specific compression analysis (Figure 11b) identifies the first soft-soil layer as the dominant contributor, accounting for approximately 55%–60% of total compression during 2025–2035. Each percentage is calculated as the modeled compression of the corresponding layer divided by the sum of modeled compression across all layers within the representative area for the same year.

The equivalent M-C formulation does not explicitly resolve viscoplastic creep, microstructural evolution, or a constitutive secondary-consolidation stage. Accordingly, the temporal and depth-dependent changes obtained from the simulations are interpreted here only as changes in modeled layer compression under the prescribed loading, groundwater, drainage, and constitutive assumptions. The contribution of the second soft-soil layer increases from 8.43% in 2025 to 22.89% in 2035, indicating progressively greater participation of deeper strata in the total modeled compression. This pattern is compatible with long-term delayed deformation of an underconsolidated soft-soil system, but it should not be interpreted as direct evidence of downward migration of a secondary-consolidation front or as confirmation of a specific creep or microstructural mechanism.

3.5.2 Anthropogenic factor: intensification of land-use intensity

While soft-soil architecture determines the baseline susceptibility to settlement, LUI introduces spatially variable surface loading. Localized overlap between higher-LUI areas and the southern subsidence zone suggests that additional development-related loading may intensify compression where it is superimposed on susceptible soft-soil deposits. Under the prescribed baseline scenario, annual LUI growth corresponds to increasing equivalent surface load and may increase modeled effective stress locally in sectors undergoing continued development. At the regional scale, the persistent south–north contrast remains primarily associated with the distribution, thickness, and consolidation state of the marine soft-soil layers, with LUI contributing to local variations in settlement magnitude.

4 Discussion

The predicted subsidence pattern indicates that mitigation in the southern high-risk sector should focus on the long-term compression of deep marine soft soils rather than on shallow foundation treatment alone. At the representative southern subsidence center, the first soft-soil layer contributes approximately 55%–60% of total compression. Where further development is planned, staged loading, surcharge or vacuum preloading with vertical drainage where applicable, and control of rapid load accumulation should therefore be considered to reduce post-construction deformation.

The baseline-scenario differential subsidence estimate exceeding 400 mm may have different implications for major infrastructure. For terminals and storage yards, uneven deformation may affect pavement serviceability, equipment alignment, and structural connections; phased construction and long-term deformation monitoring are therefore necessary. Buried pipelines may require flexible joints, settlement-compatible supports, and regular inspection at transition zones. For flood-control and drainage facilities, differential settlement may reduce effective crest elevation and drainage gradients, and should be incorporated into hydraulic-capacity assessment and maintenance planning.

At the regional scale, land-use planning should combine geological susceptibility with allowable development loading. Higher-intensity development should be carefully controlled in areas underlain by thick underconsolidated soft soils, while long-term monitoring should integrate surface settlement, layered deformation, groundwater level, and pore-water pressure. These observations can support periodic model updating and the timely adjustment of engineering and operational measures.

The modeled concentration of larger subsidence in reclaimed and intensively developed sectors is broadly consistent with InSAR observations from global coastal cities, the southwestern region of Bangladesh, and the Pearl River Delta, where spatial variability has been associated with combinations of compressible deposits, reclamation, groundwater conditions, and human activity (Wu et al., 2022; Qiao et al., 2024; Ren et al., 2024; Sarker et al., 2024). Unlike monitoring- or attribution-oriented studies, the present framework translates LUI into an equivalent mechanical boundary within a borehole-constrained three-dimensional model. This distinction supports scenario screening but does not establish a causal decomposition of geological and anthropogenic contributions.

Although the three-dimensional model retained the spatial variation in the elevations and thicknesses of the generalized stratigraphic units, uncertainties remain due to sparse borehole control, spatial interpolation, representative material parameters, the simplified groundwater boundary, and the equivalent parameterization of delayed compression. The modified M-C correction is intended for underconsolidated marine soft soils with estimated consolidation degrees of approximately 70%–80%, under geological, groundwater, loading, and stratigraphic conditions comparable to those represented in the Meishan model. Because the correction factors were applied and evaluated within the borehole-derived three-dimensional stratigraphic geometry, the correction is not formulated as a universal numerical threshold for soft-soil thickness. Recalibration is required where soft-soil thickness, consolidation state, groundwater conditions, or loading intensity differ substantially from those of the study area. These factors may affect the predicted magnitude and spatial boundary of localized settlement. The present analysis is structured as a deterministic baseline-scenario assessment rather than a probabilistic forecast. Accordingly, the 2025-2035 annual and cumulative values are conditional baseline-scenario estimates rather than confidence-bounded forecasts; plausible changes in stratigraphy, material properties, groundwater, permeability, equivalent surface loading, or development trajectory may change their magnitude and timing. The use of spatially uniform correction factors further limits the representation of local heterogeneity. The monitoring and PS-InSAR comparisons provide local and regional corroboration under the Meishan conditions but do not separately quantify the errors arising from the individual simplifications. The results should therefore be interpreted as regional screening estimates.

Future work should incorporate additional boreholes, alternative stratigraphic-discretization schemes, borehole-density and soft-soil-thickness sensitivity analyses, long-term groundwater monitoring, and calibrated M-C–SSC comparisons to quantify these uncertainties. HMG06 was used as the designated co-located multi-level settlement site for depth-specific evaluation. The remaining boreholes constrain the three-dimensional stratigraphic geometry and baseline properties, while the regional surface response is assessed through the 2022 PS-InSAR pattern. Because the LUI-to-load relationship was informed by the processed 2022 PS-InSAR product, the 2022 PS-InSAR comparison is treated as a regional spatial consistency assessment rather than fully independent validation. It uses the same correction rules across the regional model, but its surface-integrated signal does not replace evaluation at a second depth-instrumented borehole. Monitoring-informed Class C prediction studies indicate that combining depth-specific settlement and pore-water-pressure observations can strengthen parameter updating and long-term settlement prediction (Huang S. et al., 2025). Consequently, the current evidence supports regional pattern screening under the Meishan conditions, not universal coefficient transferability; a multi-site calibration-validation design with depth-resolved settlement and, where possible, pore-water-pressure observations is required before applying the numerical factors elsewhere.

5 Conclusion

This study investigated the spatiotemporal evolution and combined geological-anthropogenic controls of land subsidence in a highly urbanized coastal reclamation environment. By integrating a refined 3D geomechanical model, LUI-based surface-load inversion, underconsolidated soil parameter correction, and PS-InSAR comparison, the long-term deformation trajectory of the Meishan Bonded Port Area from 2025 to 2035 was analyzed. The principal conclusions are as follows:

• The coupled multi-factor framework links urban spatial expansion with the modeled response of underconsolidated marine soft soils. At HMG06, the annual total-compression relative errors for 2018-2021 were 2.20%–7.80%; across the 12 repeated layer-year observations, MAE was 1.22 mm, RMSE was 1.62 mm, and mean bias was −0.23 mm. These descriptive results support local plausibility at one monitored column but do not constitute independent full-domain validation.
• The southern reclamation zone remains the primary subsidence center. Under the prescribed baseline scenario, the maximum annual subsidence rate decreases from 35 to 40 mm/year in 2025 to 15–20 mm/year in 2035, a modeled trend compatible with attenuation of delayed compression rather than a direct constitutive simulation of creep. Nevertheless, cumulative deformation continues to increase, with peak cumulative subsidence modeled to reach 400–450 mm by 2035.
• Stratigraphic analysis at the representative southern subsidence center shows that the underconsolidated first soft-soil layer contributes approximately 55%–60% of total compression, indicating its modeled dominant contribution in the high-subsidence sector.
• Anthropogenic LUI intensification modifies local settlement through additional surface loading. Localized overlap between higher-LUI areas and compressible soft soils may enhance deformation, while the broader south–north pattern is mainly governed by the distribution and consolidation state of the marine soft-soil deposits. These results indicate an association within the prescribed model structure, not a quantitative causal separation of anthropogenic and geological contributions.
• The modified M-C scheme provides a site-conditioned regional screening approximation of delayed soft-soil compression when the objective is consistent regional screening rather than layer-specific creep calibration. The transferable contribution is the transparent workflow that combines borehole-constrained geometry, declared equivalent factors, depth-resolved local evaluation, and a regional spatial consistency assessment; the 20%–30% correction factors are not universal and require local recalibration and independent validation before application outside comparable Meishan conditions. The scheme is not a substitute for SSC-based site-specific design.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative artificial intelligence was used only for manuscript formatting and language-editing support. Specifically, OpenAI ChatGPT (GPT-5.5 Thinking, OpenAI) was used to assist with formatting according to Frontiers guidelines, dehyphenation, and preparation of the scope statement. The authors reviewed and verified all scientific content, references, data interpretation, tables, and figures.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Zhang Cong, Central South University Forestry and Technology, China

Reviewed by: Navakanesh M Batmanathan, National University of Malaysia, Malaysia

Tao Peng, Guizhou Communication Polytechic University, China

Md. Shahoriar Sarker, University of Dhaka, Bangladesh

Footnotes

InSAR, Interferometric Synthetic Aperture Radar; PS-InSAR, Persistent Scatterer Interferometric Synthetic Aperture Radar; LUI, Land-Use Intensity; M-C, Mohr-Coulomb; SSC, Soft Soil Creep; MAE, Mean Absolute Error; RMSE, Root Mean Square Error.