Mechanism of vertical displacement beneath shallow compression zone in coastal area

Previous studies showed that there is a discrepancy between total subsidence (measured on the ground surface) and shallow compression (measured within several hundred meters beneath the ground surface) in Taiwan. This dif‑ ference is referred to as deep displacement in this study. The variations of deep displacement are opposite to those of ground surface displacement and groundwater level within the depth of several hundred meters. The mechanism is unknown and requires further investigation. This study adopts two kinds of geodetic observation data and sepa‑ rates land subsidence into shallow compression and deep displacement to investigate the mechanism of deep displacement. A tectonically active coastal area in Taiwan is selected as the study area. The assessment results show that the associated variations are likely due to cyclic hydraulic loading and unloading. The variations of deep dis‑ placement are opposite to those of ground surface displacement. This study proposes that these variations are due to hydraulic expansion and contraction. The mechanism is demonstrated using a hydromechanical model. The results of a cross‑correlation analysis show that hydraulic expansion and contraction occur at certain depths. The study results provide important information on the mechanism of deep displacement that can be used in tectonophysical and land subsidence investigations.


Introduction
Mechanisms induced by natural factors, such as natural compression, seasonal hydrological loading, tectonic activity, fault activity, and earthquakes, contribute to land subsidence.Fu et al. (2012), Chanard et al. (2014), and Hsu et al. (2020) reported that the seasonal hydrological loading of continental water storage and its loading effects are significant in southern Alaska, the Himalayas, and Taiwan, respectively.Hsu et al. (2021) found that the hydrological cycle in western Taiwan causes elastic loading and unloading effects, which may be the primary driving mechanism of the observed synchronized modulation of earthquakes.In the polar region or highelevation areas in winter, frozen water and snow create mass loading on the ground surface, which can lead to land subsidence.The recovery of frozen soils in summer mainly depends on their physical properties (Martens 2016;Knappe et al. 2019).These mechanisms mainly induce land subsidence at shallow depths (on the scale of hundreds of meters).Carminati and Di Donato (1999) and Teatini et al. (2007) reported that there is a small natural subsidence in Italy (− 1.6 and − 1.2 mm/yr in Ravenna and Venice, respectively).Land subsidence can also be induced by tectonic activity, fault activity, creep, or continuous plate and fault movement (Wu 1978;Yu et al. 1997;Lacombe et al. 2001;Bos et al. 2003;Chiang et al. 2004;Shyu et al. 2005;Hu et al. 2006Hu et al. , 2007)).Parsons (2006) reported that earthquakes result from elastic stress built up in part by tectonic motion.Earthquake events can induce a sudden change in displacement and tectonic activity and fault activity can induce a long-term movement of a plate or geological structure, inducing strata displacement in both the vertical and horizontal directions.These mechanisms mainly induce land subsidence at deep depths (on the scale of kilometers).However, the assessments in these studies did not always consider the influence of anthropogenic activity, especially in areas with dense populations, which resulted in large uncertainty in the assessment results.
In some delta regions, such as the Mekong Delta in Vietnam (Zoccarato et al. 2018), the Po Delta in Italy (Teatini et al. 2011), and the Mississippi Delta in the USA (Törnqvist et al. 2008;Jankowski et al. 2017), young sediments (e.g., Holocene sediments) commonly have a high compression rate because of a relatively loose sedimentary environment and high groundwater over-extraction.Liu et al. (2020) assumed that land subsidence has two components: (1) a contribution from bedrock systems or non-compressing strata caused by tectonic subsidence and (2) a contribution from the compression of compressible, non-bedrock aquifer systems caused by primary consolidation due to subsurface fluid withdrawal.Therefore, coastal areas with anthropogenic activity might have serious land subsidence due to shallow compression induced by construction loading and/or groundwater over-extraction (Poland 1984;Galloway et al. 1998;Alley et al. 2002;Hou et al. 2005;Gambolati et al. 2006;Hu et al. 2006;Meckel et al. 2006;Hsu et al. 2007;Teatini et al. 2007;Hung et al. 2010;Galloway and Burbey 2011;Hsieh et al. 2011;Tung and Hu 2012;Erban et al. 2013;Wang et al. 2015aWang et al. , 2015b;;Hwang et al. 2016;Jones et al. 2016;Minderhoud et al. 2017;Tran et al. 2022).The depth at which anthropogenic activity-induced compression has influence is relatively shallow, but the contribution to land subsidence is relatively large compared to that of tectonic activity and plate displacement.Tectonophysical studies commonly use ground surface observation systems to assess tectonic activity and plate displacement (Yu et al. 1997;Bos et al. 2003;Hou et al. 2005;Shyu et al. 2005;Hu et al. 2006Hu et al. , 2007;;Parsons 2006;Ching et al. 2007Ching et al. , 2011;;Wang and Shen 2020).Therefore, the influence of anthropogenic activity can have a large effect on the investigation of mechanisms in tectonophysical studies.
Techniques such as leveling, the Global Positioning System (GPS), deep benchmarks, multi-layer compaction monitoring wells (MLCWs), and interferometric synthetic aperture radar (InSAR) are commonly used to monitor land subsidence.Leveling, GPS, and InSAR are used to monitor total subsidence on the ground surface and deep benchmarks and MLCWs are used to monitor compression within a certain depth range.Previous studies (Hung et al. 2010(Hung et al. , 2018) ) showed that there is a discrepancy between total subsidence (measured by leveling and GPS) and shallow compression (measured by MLCWs) in Taiwan.This difference is caused by displacement beyond the major compression zone, which is referred to as deep displacement in this study.The mechanism of deep displacement is unknown due to a lack of direct observation data.Deep compression due to groundwater over-extraction has been proposed as a possible mechanism (e.g., Galloway et al. 1998;Erban et al. 2013;Hung et al. 2018).However, the variations of deep displacement are opposite to those of groundwater level variations at shallow depths (Tran and Wang 2020), which is inconsistent with deep compression due to groundwater over-extraction (Lu et al. 2020;Hsu et al. 2021).Further investigations are thus needed.
The coastal area of Pingtung Plain in Taiwan had the most substantial cumulative subsidence (3.56 m), which was mainly due to groundwater over-extraction, in the period 1972-2018 (Taiwan Water Resources Agency (WRA) 2018).Therefore, the WRA installed several subsidence and groundwater monitoring systems in this area.Previous studies reported that approximately three-quarters of land subsidence in Pingtung Plain is caused by groundwater over-extraction and that part of the remaining land subsidence is caused by tectonic activity (Kuo et al. 2001;Hou et al. 2005;Hu et al. 2006Hu et al. , 2007)).However, these studies did not provide observation data to quantify land subsidence due to anthropogenic activity and natural factors and did not discuss the subsidence mechanisms at shallow and deep depths.Therefore, Pingtung Plain is a suitable location for validating the mechanism of vertical displacement due to natural (deep depth) and anthropogenic (shallow depth) factors.The present study uses multiple geodetic techniques to quantify and investigate land subsidence due to deep displacement in a tectonically active region.A coupled hydromechanical model is used to demonstrate the phenomenon of hydraulic expansion and contraction, defined in this study, for deep displacement.

Study area
The study area is located in the coastal area of Pingtung Plain, Taiwan, for which there are sufficient data to investigate the mechanism of land subsidence (Fig. 1).A borehole log and a hydrogeological cross section between three wells, namely Dexing, Fangliao, and Jiadong wells, in the study area were constructed, as shown in Fig. 2A and B. According to the geological characteristics of the alluvial fan, the study area stretches from the proximal fan (Fangliao area) to the distal fan (Jiadong area).The energy available for the sediment process is low due to the apparent terrain elevation from the northeast to the southwest of Linbian River (Fig. 1B).The geological age has a marked change from 11,030 ± 60 to 29,200 ± 330 years ago; thus, the Holocene unconsolidated sediment is distributed within a depth of about 120 m (Fig. 2B).S-wave velocity and electrical resistivity data provide information on the lithology at Jiadong well, as shown in Fig. 2C.These data can be used as a reference for determining the depth of cemented sediment.
The tectonic processes in Pingtung Plain have two states (Fig. 1A), namely a dynamically extruded state in the east and southeast and a static state in the plain area, separating Pingtung Plain by the Chaochou active fault (Lacombe et al. 2001).There is a tectonic escape connected with a tectonic extrusion and a lateral extrusion that operates at the southern tip of the Taiwan collision belt, towards the submerged area in the Manila trench, and acts as a free boundary (Lacombe et al. 2001;Hou et al. 2005;Hu et al. 2006Hu et al. , 2007;;Ching et al. 2007).Quaternary unconsolidated sediments form a basin and lie on the wedge-top depozone and associated piggyback basins (Chiang et al. 2004).They are caused by the weathering of the Taiwan orogeny and form the principal aquifer of Pingtung Plain.(1997) and Ching et al. (2007).B Distribution of land subsidence monitoring stations and groundwater level monitoring wells, where fault system is from Chiang (1971) and Ching et al. (2007).

Monitoring systems and observation data
In Taiwan, four techniques are commonly used to observe land subsidence, namely leveling, GPS, MLCWs, and InSAR (Hung et al. 2010(Hung et al. , 2021)).MLCWs, developed in Taiwan, are used to measure the compression of the aquifer system in different layers (Hung et al. 2010).In this study, Jiadong and Fangliao MLCWs are adopted due to their long monitoring period and proximity to GPS stations, which measure deformation within 200 m from Earth's surface.GPS data are adopted due to their high sampling rate, which allows a displacement comparison at a given time point between two observation systems.Two GPS stations, named CLON and FALI, are selected in this study because they are the closest stations to the MLCWs.Furthermore, five monitoring wells, named Dazhuang (DZ), Wenfeng (WF), Daxiang (DX), Fangliao (FL), and Dexing (DeX), respectively, are used for groundwater level observation in the study area.The distribution of the observation systems is shown in Fig. 1B.Note that the observation stations are not completely collocated.It is assumed that for a given station, the closest station records similar behavior.In this study, the behaviors recorded at nearby stations are compared.The comparison of noncollocated data is a source of uncertainty in this study.
The observation data of hydrology and land subsidence in the study area from 2007 to 2016 are shown in Fig. 3.The variations of land subsidence, groundwater level, and rainfall quantity show similar patterns, which are induced by seasonal hydrological variations.The dry season is from November to April, during which there is a steady decline in the groundwater level due to small recharge and large discharge.Both anthropogenic (e.g., pumping for drinking water, irrigation, and aquaculture) and natural (e.g., evapotranspiration, submarine groundwater discharge, and baseflow) factors decrease the groundwater level, inducing the compression of porous media.The wet season is from May to October, during which precipitation is mainly from monsoon rains, plum rains, and typhoons, which largely increase groundwater recharge.Hence, groundwater level increases, leading to the dilation of porous media.Figure 4 shows the seasonal fluctuations of groundwater level in three aquifers (defined in Fig. 2B).The fluctuations are 0.82 − 32.78, 11.79 − 33.47, and 7.44 − 22.63 m for the first, second, and third aquifers, respectively.A good permeable aquifer lets water easily recharge into and discharge from the aquifer and induces a large variation of groundwater level.

Deep displacement calculation
Based on the monitoring limit of MLCWs, vertical displacement was separated into a shallow component (within 200 m) and a deep component (beyond 200 m).In this study, the former is denoted as shallow compression and the latter is denoted as deep displacement.Note that this study does not attempt to define the threshold of deep displacement; it only collects observation data to investigate the mechanism of vertical displacement beneath the major compression zone.Different areas may have different thresholds for the vertical displacement, with different responses above and below the threshold.
To obtain the accurate results of deep displacement, both the GPS and MLCW data need to be assessed carefully.Satellite signals (e.g., GPS) are adversely influenced by the atmosphere and air pressure variations.It is difficult to avoid numerical errors when using real-time resampling.The measurement period of MLCWs is about 1 month, whereas that of continuous GPS is 1 day.Therefore, the weekly average from the mid-point solution of GPS data is calculated to reduce numerical errors.The measurement date of MLCW data is the reference point.GPS data collected within three days before and after the reference date of MLCW data are averaged (for a total of 7 GPS data points).The MLCW data obtained at the same time point are subtracted from the averaged results of GPS data to calculate their difference (i.e., the deep displacement).Deep displacement is calculated at each time point of MLCW data to obtain the time series of deep displacement.

Hydromechanical model
A coupled hydromechanical model is used to demonstrate the phenomenon of hydraulic expansion and contraction proposed in this study.The commercial software COM-SOL Multiphysics and poroelastic theory are adopted for the simulation.The governing equation for fluid mass conservation can be written as: where ρ f is the fluid density, S is the storage model (writ- ten as , where c b is the com- pressibility of the bulk material, α is Biot's effective stress coefficient, c f is the compressibility of the fluid, and n is the porosity), p is the change in pore water pressure, k p is the permeability, µ is the fluid dynamic viscosity, ε v is the volumetric strain, and Q f is the sink/source term.
The governing equation for force equilibrium of solid mechanics and the stress-strain constitutive relationship with linear material properties can be written as: (1) (2) where σ ij is the Cauchy stress tensor, I is the identity matrix, C ijkl is the elastic matrix ( , where E is Young's modu- lus, υ is Poisson's ratio, and , and ε kl is the Cauchy strain tensor.The strain ε ij can be expressed as: where ∇u is the displacement gradient tensor and T denotes the transpose of the matrix.
To demonstrate the mechanism of hydraulic expansion and contraction, a simplified hydrogeological profile based on the Fangliao MLCW was developed, as shown in The setting parameters and inverted results are listed in Table 1.Note that not all the values of these parameters are from experiments; some were obtained from the literature and model inversion because most of the required parameters are unavailable, especially the properties of the aquitard.
The initial condition is assumed to be the steady-state condition with zero change in pore water pressure and no displacement in the whole domain.Thus, a positive (negative) change in hydraulic head represents the hydraulic head being higher (lower) than the initial value.The sink and source due to seasonal variations are the driving force in the numerical model.Twelve months with seasonal variation, namely January to April (dry season), May to October (wet season), and November and December (dry season), based on the hydrological conditions of the DeX groundwater observation were used to simulate the hydraulic expansion and contraction.In the dry (wet) season, the hydraulic head decreases (increases) due to groundwater discharge (recharge).Aquifer 1 is an unconfined aquifer and thus the change in hydraulic head is smaller than those for aquifers 2 and 3, which are confined aquifers.

Data analysis
The rainfall quantities are not uniformly distributed in southern Taiwan.Annual rainfall is approximately 2380 mm, 95% of which is in the wet season.Groundwater level in the study area varies largely between the wet and dry seasons and induces an obvious variation of surface displacement (Figs. 3 and 4).In the wet season, the annual expansion of unconsolidated sediment at a shallow depth is 19.1 and 11.3 mm in the Jiadong and Fangliao areas, respectively, and the annual compression in the dry season is − 24.5 and − 16.5 mm in these areas, respectively.The compression over time in the Fangliao area showed a gradual decrease in subsidence rate after 2013, but the subsidence rate remained stable in the Jiadong area.The difference between the two areas may be attributed to differences in hydrogeological material deposition in different regions of the alluvial fan and differences in the compression stage of soil, including the time-delay phenomenon of clay compression.Land subsidence occurs in the dry season but rebounds in the wet season, which is consistent with the rainfall quantities.The rebound of shallow strata in the wet season is less than the compression in the dry season, leading to a longterm subsidence trend, which has been reported to be mainly due to groundwater over-extraction (Hung et al. 2010(Hung et al. , 2018)).The long-term shallow compression rates calculated from the MLCW data are − 9.7 and − 6.7 mm/ yr in the Jiadong and Fangliao areas, respectively.

Correlation analyses
The calculated correlation coefficients between observation data are listed in Table 2. Time lag was not obvious in the correlation calculations and is thus not shown here.The correlation coefficients between the MLCW data and groundwater levels at different wells are high.Because this area is subjected to serious groundwater extraction, the high correlation implies that the shallow compression is mainly due to groundwater over-extraction induced by anthropogenic activity.The correlation coefficients between groundwater levels at different depths at a given location are high, which indicates that the variation patterns at different depths are consistent due to groundwater discharge in the dry season and rainfall recharge in the wet season.Fourier analysis shows that all these data have predominant periods of about 1 year.Therefore, the rainfall, shallow groundwater level, deep groundwater level, and MLCW observations indicate that the driving forces of the variations in these time series are the same.However, after the removal of shallow compression, the correlation coefficients between deep displacement and other observations become very small and even negative, which implies that the mechanism of deep displacement is different from that of shallow compression.

Deep displacement
Figure 6A and B shows the calculation results of deep displacement in the Jiadong and Fangliao areas, respectively, along with shallow compression and rainfall quantities.The contribution of deep displacement to total land subsidence is 37.06% and 13.63% in the Jiadong and Fangliao areas, respectively (Fig. 7).The average rates of deep displacement are − 4.6 and 0.9 mm/yr in the Jiadong and Fangliao areas, respectively.Deep displacement in both areas shows fluctuation over time.The predominant periods are about 1 year, which might be related to the hydrological cycle.However, most of the correlation coefficients between groundwater level and deep displacement are low or even negative (Table 2).The calculated correlation coefficients between deep displacement and rainfall (calculated as cumulative rainfall departure; Weber and Stewart 2004) in the Fangliao and Jiadong areas are both negative (− 0.47 and − 0.10, respectively) and exhibit no time lag.Although the correlation of shallow compression between the Jiadong and Fangliao areas is high (0.95), the trends of deep displacement in these two areas are totally different; their correlation coefficient is only 0.42 (Table 2).
After the removal of shallow compression, the correlation coefficients between deep displacement and other observations become very small and even negative (Table 2), which implies that the mechanism of deep displacement is different from that of shallow compression.To understand the mechanism of deep displacement, the threshold between shallow and deep components was set to various depths.For each station, three depths were set as the threshold to calculate shallow compression and deep displacement.The MLCW data for depths of 47, 118, and 194 m in the Jiadong area and those for depths of 67, 121, and 191 m in the Fangliao area were selected based on the MLCW measurement depth and the hydrogeological model of the aquifer system in Pingtung Plain.For each depth, shallow compression was measured based on MLCW data from the land surface to the given depth and deep displacement was calculated based on the difference between GPS values and MLCW values.The results of cross-correlation are listed in Table 3.The time lag information calculated from cross-correlation is also shown in the table to demonstrate that the time lag is not obvious.It is interesting to see that the deep displacement at shallow depths has high correlations with hydrological observations but the quantities decrease with increasing depth.The results clearly show that at a certain depth, the correlation between deep displacement and hydrological observations transitions from positive to negative.

Discussion
The deep displacement (beyond a depth of 200 m) results in both areas fluctuate temporally, which might be related to the hydrological cycle.The mechanism of the variations of deep displacement could be cyclic hydraulic loading (e.g., Fu et al. 2012;Chanard et al. 2014;Martens 2016;Knappe et al. 2019;Hsu et al. 2020Hsu et al. , 2021;;Tran 2020), which induces the hydraulic expansion and contraction proposed in the present study (Fig. 8).Cyclic hydraulic loading is an additional form of seasonal loading and unloading and contributes to the total mass variations at a certain depth.An increase (decrease) in pore water pressure in the aquifer induces an expansion (contraction) condition that both uplifts (sinks) the top strata Table 2 Correlation coefficients for various components JD-DD: deep displacement in Jiadong area; FL-DD: deep displacement in Fangliao area; CRD: cumulative rainfall departure.The abbreviations of monitoring wells are defined in Figs. 2 and 3. Note that this table is separated into two parts.The first three columns and rows assess the correlation between deformation in different areas with hydrological data; thus, the quantities are not symmetrical between columns and rows.The fourth to fourteenth columns and rows show the correlation between the hydrological data; thus, the quantities are symmetrical between columns and rows  and loads (unloads) the bottom strata.The effect of an increase in pore water pressure on the top strata is an uplift of the ground surface and that on the bottom strata is hydraulic loading.The groundwater level in the wet season at a deep aquifer (e.g., DZ2, DeX3, or FL2) can increase pore water pressure at a depth of 200 m (Fig. 4C), which is large enough to create observable hydraulic loading on strata deeper than 200 m.A negative (positive) value of annual displacement indicates the downward (upward) movement of deep strata due to increasing (decreasing) pore water pressure.Figure 9 shows the phenomenon of cyclic hydraulic loading and unloading.In Fig. 9A, detrended deep displacement and detrended hydraulic factors show an obviously negative relationship in the Fangliao area.In Fig. 9B, during the wet season, the average downward displacements are − 1.6 and − 14.6 mm/yr in the Jiadong and Fangliao areas, respectively, and during the dry season, the average uplift displacements are 5.9 and 12.6 mm/yr in these areas, respectively.Cyclic hydraulic loading and unloading are obvious in the Fangliao area but not in the Jiadong area,   A coupled hydromechanical model is used to demonstrate the phenomenon of hydraulic expansion and contraction.The simulation results are shown in Fig. 10.The phenomenon of hydraulic expansion and contraction can be demonstrated by the simulated strain in the developed hydromechanical model.The simulated strain in aquifers 1 to 3 shows a decreasing trend (compression) in the dry season and an increasing trend (dilation) in the wet season, which is consistent with observations.However, the simulated strain in aquitard 3 shows the opposite responses, that is, an increasing trend (dilation) in the dry season and a decreasing trend (compression) in the wet season.This reversed phenomenon is the cyclic hydraulic expansion and contraction in the aquifer system proposed in this study.Note that this phenomenon only occurs when the formation has a small permeability (Table 1).Such conditions include, for example, the presence of an aquitard capable of forming an interface to absorb stress caused by variations in pore water pressure within an adjacent aquifer.Consequently, the occurrence of hydraulic expansion and contraction is contingent upon the strata meeting specific criteria, such as possessing a very small permeability at certain depths.The findings from the hydromechanical simulations corroborate observed behaviors and support the mechanism proposed in this study.
The linear trend of deep displacement in Fig. 6 might be due to tectonic activity, as discussed in the literature (e.g., Kuo et al. 2001;Hou et al. 2005;Hu et al. 2006Hu et al. , 2007)).Figure 11 shows that the direction of tectonic escape is toward the southwest, with an azimuth that ranges from 243.9° to 245.7° (Tran and Wang 2020).In this study, an average horizontal displacement rate of 26.6 mm/yr was calculated using GPS data; this value is consistent with  Lacombe et al. 2001;Hou et al. 2005;Hu et al. 2006Hu et al. , 2007;;Ching et al. 2007).The observation stations correspond to those in Fig. 1 the findings of Hu et al. (2006).Furthermore, Hu et al. (2007) found that a tectonic escape connected with a tectonic extrusion and a lateral extrusion operates at the southern tip of the Taiwan collision belt.It is thus important to recognize that the horizontal displacement of a stratum corresponds to its vertical displacement since the extension and contraction of the stratum can also change its thickness or elevation (Hu et al. 2007).However, since there is no further evidence to support the mechanism of the linear trend of deep displacement, it is not discussed in this study.

Conclusions
GPS, MLCW, rainfall, and groundwater level data were collected to investigate the mechanism of deep vertical displacement in Pingtung Plain, Taiwan, which is located in a tectonically active coastal area.The collected data all showed obvious seasonal variations with a predominant period of about 1 year.The variations of deep displacement also have a 1-year period but are opposite to those of hydrological data, indicating that these variations are likely due to cyclic hydraulic loading and unloading.The large variations of groundwater level at a depth of 200 m create cyclic hydraulic loading and unloading for strata at depths deeper than 200 m, which explains the variations of deep displacement.The correlation coefficients between hydrological data and deep displacement at various depths show a decrease with depth (without obvious time lag), from a high positive correlation to a negative correlation.The results imply that the phenomenon of hydraulic expansion and contraction can be observed at certain depths.A simulation using a coupled hydromechanical model was conducted.The simulation results demonstrate that the phenomenon of hydraulic expansion and contraction may occur when a deep formation has a small permeability.The numerical model results are consistent with the observations and the proposed mechanism.In the Jiadong and Fangliao areas, deep displacement contributes 37.06% and 13.63% to total subsidence, respectively.These contributions are similar to that for tectonic activity estimated in the literature for the study area, indicating that the linear trend of deep displacement might be due to tectonic activity.However, there is no evidence to support this mechanism; further investigations are thus needed.
The variations of deep vertical displacement are opposite to those of ground surface vertical displacement.The opposite responses indicate that the effects of natural factors on the ground surface obtained without the removal of the influence of anthropogenic activity will lead to large uncertainty in the assessment results.The study results provide important information on the mechanism of deep vertical displacement, which support the suggestion, proposed in the literature, that groundwater level variations might trigger earthquakes.This study provides partial evidence for the mechanism of deep vertical displacement.Further observations and investigations are required to provide more evidence.

Fig. 2
Fig. 2 Geological information.A Borehole log, B cross section of hydrogeology, and C S-wave velocity (Kuo et al. 2016; Chen et al. 2022) and electrical resistivity geophysical survey at Jiadong area with lithology from Geological Survey and Mining Management Agency (GSMMA).The trace of the hydrogeological profile and S-wave station are shown in Fig. 1

Fig. 3 Fig. 4
Fig. 3 Observation data for vertical displacement and groundwater levels in A Jiadong and B Fangliao areas.The MLCW data in the period 2013-2014 are missing.A negative (positive) value indicates the compression (expansion) of unconsolidated sediment Fig. 5.According to the geological setting in the Fangliao area, three aquifers and three aquitards were set within a depth of 200 m.The bottom layer (aquitard 3) is considered to demonstrate deep displacement caused by hydraulic loading and unloading.The central column, with a width of 100 m, was the main area used to simulate hydraulic expansion and contraction.The boundary of the two sides of the hydrogeological profile was set to have no horizontal deformation and zero change in pore (3) ε ij = 1 2 ∇u + (∇u) T , water pressure.To mitigate the boundary effect on the main area, the modeling domain was extended to 2.5 km away from the central column.The bottom boundary was set as a no-flow boundary with zero (fixed) displacement.The top boundary was set to have zero change in pore water pressure and the free-traction condition.A negative (positive) displacement indicates a downward (upward) movement or a compression (dilation) of strata.Six observation points in the central column (marked in Fig. 5) were chosen to demonstrate the simulation results.

Fig. 6
Fig. 6 Comparison between observed subsidence and calculated deep displacement with rainfall in A Jiadong area and B Fangliao area

JDFig. 8 Fig. 9 A
Fig. 8 Conceptual model of deep displacement induced by hydraulic loading and unloading.Downward (upward) movement of deep strata is related to hydraulic loading (unloading) during wet (dry) season

Fig. 10
Fig. 10 Simulated A change in pore water pressure in three aquifers and B strain in three aquifers and aquitard 3 (i.e., bottom layer)

Table 1
Material properties for hydromechanical model

Table 3
Correlation coefficients for various deep displacement depths