Model performance
The performance of SR such as the RMSEs in the multi- or single-aquifer
system (RMSE = 0.6–0.8 cm) is better than that of OLS (RMSE =
1.58–2.11 cm) (Table 1). The best models are appear to be Model 1 in
the multi-aquifer system and Model 4 in the single-aquifer system.
However, the subsidence–drawdown model using OLS is not appropriate due
to its poor performance (Table 1). The SR is more accurate than that of
OLS. Results of OLS reveal a subsidence–drawdown function with negative
a slope for Layer 1 (high drawdown -low subsidence exists in the
proximal fan) but a positive slope for Layers 2, 3, and 4 in Model 1
(Table 2). Compared with Models 2–5, the best SR with the lowest RMSEs
occurs in Model 4 because the subsidence–drawdown correlation is the
highest, but the poorest model with highest RMSE is Model 2 due to its
low subsidence–drawdown correlation (Fig. 3). Considering the
validation using the GPS data of ground surface (Fig 1), the RMSEs of
Model 4 using linear regression and SR are 1.52 and 0.67 cm,
respectively. Validation of our subsidence model with available GPS data
provides consistent and accurate results.
Fig. 5 shows the total estimated land subsidence calculated from linear
regression and SR using all-layers (Model 1) and single layer drawdowns
(Models 2–5). In a linear regression, the estimated land subsidence
bowl occurs in the distal fan area (near the coast), but also in
the proximal
and distal fan areas in Model 2 due to the variable drawdown patterns
occurring in Layer 1. The estimated maximum subsidence in the OLS models
occurs at the boundary area of Yunlin, whereas in the SR models
subsidence occurs in the inland area of Yunlin (Fig. 5). However, the
estimated pattern using OLS does not match the observed subsidence. The
SR explicitly considers spatially dependent models to overcome the
spatial variability of land subsidence mapping because the SR extends
the linear regression by estimating a set of spatial parameters rather
than one single set of parameters (Brunsdon et al., 1996; Fotheringham
et al., 2003). The model RMSEs are less than or equal to 0.8 cm so that
the spatial regression can be used to model spatially varying
relationships between subsidence and drawdown in the aquifer system.
Spatial model coefficient
and implication
Fig. 6 exhibits the spatial pattern of regression slope coefficients of
the SR in Layers 1–4 for Models 2–5. The regression coefficient
distribution varies spatially. The patterns of spatial coefficients are
typically similar in Layers 2, 3, and 4. Most coefficients have a
positive slope for total drawdown and subsidence (Fig. 6). Only the
proximal fan and the area close to Wu River exhibits a negative relation
between the total drawdown and land subsidence (Fig. 6). A negative
relation can occur where unloading causes a small soil rebound or where
drawdown continues but the subsidence rate declines due to recharge
effects (i.e., rainfall, ponds/lakes, canals, and irrigation).
The long-term trend of hydraulic head and vertical displacements are
useful to estimate the inelastic skeletal storage coefficient for
subsidence features (Miller and Shirzaei, 2015). The spatial patterns of
skeletal storage coefficients in four layers can be estimated using the
slope coefficients of spatial regression (Fig. 6), which also reflect
the skeletal storage coefficients
of the aquifers. In Fig. 6(c), the large coefficient value (purple
color) represents inelastic compaction in the inland area, whereas the
small coefficient value (blue color) represents elastic compaction in
other areas. The resulting subsidence pattern clearly suggests
elasto-plastic mechanical behavior, which includes elastic and inelastic
skeletal storage in over-consolidated and normally consolidated soils
(Hung et al., 2012; Liu et al., 2004). In recent decades, the excessive
pumping of groundwater has resulted in serious subsidence problems along
the coast. This subsidence coincides with the increase of aquaculture
fisheries since the 1970s. Currently, the rate of compaction in the
coastal area is declining. However, the subsidence bowl is moving inland
due to over-pumping in this region (Hwang et al., 2008). Interferometric
Synthetic Aperture Radar (InSAR) offers a powerful and useful method for
identifying regional land subsidence (Bürgmann et al., 2000; Bonì et
al., 2016; Chaussard et al., 2014b; Galloway et al., 1998b; Hoffmann et
al., 2001). The present findings match the results from previous InSAR
interferograms for the study region (Tung and Hu, 2012). Moreover, the
advent and use of GPS has made land subsidence monitoring more efficient
and cost effective. However, GPS is limited because it measures only the
surface deformation at point locations (Hosseini et al., 2007). The
approach developed in the present study provides robust and accurate
maps of estimated land subsidence that are far more spatially
descriptive than those that GPS typically provides.
In the developed models, spatial patterns of subsidence are estimated
solely on the basis of drawdown information. The effect of
groundwater-level changes on land subsidence can be modeled using a
stress–strain relationship with the elastic or inelastic skeletal
storage coefficient of the underlying aquifers. However, the entire
system includes leveling surveys, borehole extensometer data, and
multilayer monitoring of groundwater levels, with the aim to better
understand the hydrological and mechanical processes of the aquifer
system and to characterize the spatial distribution of land subsidence
with a traditional numerical model (Zhang et al., 2014). Characterizing
the hydrogeological setting of the aquifer system is a critical step to
accomplish this goal (Liu et al., 2004). The observational data used to
characterize hydrogeology is limited (Hubbard and Rubin, 2000) including
the hydraulic diffusivity of the confining layer, the distance from the
pumping wells, length of the recovery cycle, and the thickness of
confining layers. All of these factors can influence the shape of
stress–strain hysteresis loops (Burbey, 2001).
In this study, further calibration of the model is not required, unlike
other numerical models that typically require an inverse modeling
component for calibration purposes. The fundamental benefit of inverse
modeling is its ability to automatically estimate parameter values that
produce the best fit between observed and simulated hydraulic heads and
flows (Franssen et al., 2009; Poeter and Hill, 1997). However,
limitations include model uncertainty (hydrogeologic characterization of
the system), and measurement and parameter uncertainties (Højberg and
Refsgaard, 2005). For example, the numerical model contains parameter
uncertainty due to the unknown pumping rates, which can limit the
ability to adequately characterize key hydrologic parameters. Therefore,
this method can be used
efficiently for subsidence management and might be widely used for
subsidence prediction solely based on drawdown. Subsidence is influenced
by drawdown due to persistent pumping and effects of dry–wet cycles and
seasonal variations of groundwater levels (Kouda et al., 2015). This
subsidence–drawdown model can be relevant in terms of policy or land
subsidence regulation.
Conclusions
This study implements a spatial regression-based subsidence-mapping
scheme on the basis of groundwater-level observations and offers an
effective method to explore the spatial patterns of subsidence through
its correlation with drawdowns, without the use of complex hydrogeologic
models. The model developed here is based on a SR method and is easy to
fit between the total drawdown and subsidence observations. A
goodness-of-fit better than OLS (RMSEs are less than or equal to 0.8 cm)
is achieved. Without requiring a more detailed hydrogeologic
investigation and extensive calibration, this model offers reasonable
detail regarding the spatial patterns of land subsidence using a
non-linear and spatially varying relationship between the water-level
observations and resulting land subsidence. This model can be relevant
for water resource policy or land subsidence regulation.
The developed model accurately estimates the spatial pattern of land
subsidence. The estimated subsidence bowl occurs in the inland area of
Yunlin, similar to the actual observed subsidence bowl location.
However, the drawdown cone occurs in the coastal area west of the
subsidence bowl. The model can estimate the spatial pattern of
subsidence that is highly related to the hydrogeological properties and
drawdown history, which impacts the location of the observed subsidence
bowl. Results also show that the model coefficients accurately reflect
the skeletal storage coefficients of the aquifer. The large coefficient
value represents inelastic compaction in the inland area, whereas the
small coefficient value represents the elastic compaction that occurs
along the coastal area. The coefficient patterns represent elastic and
inelastic skeletal storage coefficients in over-consolidated and
normally consolidated soils.
In a future study, the recharge
effects and hydrodynamic lag that exists between the drawdown and
subsequent subsidence can be considered. The prediction of long-term
subsidence or subsidence time series can then be investigated.
Data Availability Statement
The data that support the findings of this study are available from
Water Resources Agency of Taiwan. Restrictions apply to the availability
of these data, which were used under license for this study. Data are
available from Chu et al. with the permission of Water Resources Agency
of Taiwan.
Acknowledgements
The authors thank the MOST for their financial support
(105-2621-M-006-011-).
Table 1. Model performance comparison for OLS and SR (Model 1: Layers
1–4; Models 2: Layer 1; Models 3: Layer 2; Models 4: Layer 3; (e)
Models 5: Layer 4)