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)