Required data for this investigation
This investigation uses monthly groundwater level observations in 2015. The water-level observations of aquifer layers 1, 2, 3, and 4 obtained from 33, 40, 28, and 28 monitoring wells, respectively, are evenly distributed over the entire Choshui River alluvial fan (Fig. 1a).
In addition, annual subsidence rates for 2015 from 662 leveling points (Fig. 1a) have been collected. A leveling network amounting to over 1,000 km in length is used to calculate subsidence for every 1.5 km interval along the leveling routes. Leveling specifications satisfy a loop closure of less than 3\(\sqrt{K}\) mm, where K is the length of the leveling circuit in kilometers. The vertical accuracy of leveling data is generally within 1 cm (Hung et al., 2010). Leveling has a high degree of accuracy but is time consuming and expensive compared with GPS. The leveling and groundwater level data are obtained from the Water Resources Agency of Taiwan.
Data Analysis Methods
Model development
Based on Terzaghi’s principle of effective-stress (Hoffmann et al., 2003a), the nonlinear stress–compaction relation is typically linearized with respect to preconsolidation stress (Shen and Xu, 2011). The calculation of total aquifer system compaction (subsidence) can be linearized where incremental changes in effective stress (hydraulic head) are typically small, as
\(\Delta z=S_{k}\Delta h\), (1)
where Δz is the subsidence, Sk is the skeletal storage coefficient, and \({h}\) is the change in hydraulic head. If the effective stress remains less than the preconsolidation stress (or less than the past maximum drawdown), the compaction is entirely elastic. In this study, the soil rebound (volume expansion after removal of mechanical stress) is disregarded from \({h}\), which represents the annual cumulative drawdown or the total annual change in hydraulic head (Chu and Chang, 2009). Rebound tends to be less than 10% of the original subsidence (Chaussard et al., 2014a). The total annual drawdown is then stated as:
\(\Delta h=\sum_{t}{h}_{t}\), (2)
where \({h}_{t}=h_{t+1}-h_{t}\) if . In this one-year case, the time step t varies from 1–12 representing the months of the year.
Subsidence is a function of accumulated drawdown. Furthermore, a regression expression identifies the relation between subsidence and groundwater level. Equation (1) is extended to yield an estimate of subsidence at observation i in a multi-layer aquifer, which can be expressed as:
\({z}_{i}=\sum_{l}{S_{\text{kl}}{h}_{\text{il}}+\varepsilon_{i}}\), (3)
where l is the layer number in the multi-layer aquifer,\(S_{\text{kl}}\) is the l th skeletal storage coefficient in a multi-layer system and \({z}_{i}\) is the land subsidence at observationi . This study uses the annual subsidence rate for 2015 and\(\ {h}_{\text{il}}\) is the estimated total groundwater decline (drawdown) at observation i from aquifer layer l .\(\varepsilon_{i}\) is the residual of the regression model. Equation (3) is further extended to allow for spatially varying subsidence and drawdown in a multi-layer aquifer system as follows:
\({z}_{i}=\sum_{l}{S_{\text{kl}}\left(u_{i},v_{i}\right){h}_{\text{il}}+\beta_{0}\left(u_{i},v_{i}\right)+\varepsilon_{i}}\), (4)
where \(S_{\text{kl}}\left(u_{i},v_{i}\right)\) varies with the spatial coordinates (\(u_{i},v_{i}\)) at observation i and is the slope of the spatial regression parameters in the l th layer.\(\beta_{0}\left(u_{i},v_{i}\right)\) is the intercept at observationi of the spatial regression parameters in the multi-layer system. The slope coefficient in Equation (4) without the intercept\(\beta_{0}\) represents the skeletal storage coefficient \(S_{k}\) in each aquifer, which is required to accurately estimate subsidence from drawdown.
Spatial coefficient estimation
The SR is an extension of ordinary least squares (OLS) (Fotheringham et al., 2003). Using the SR, the estimated parameter matrix\(\hat{\beta}\left(u_{i},v_{i}\right)\ \), which includes\(S_{\text{skl}}\left(u_{i},v_{i}\right)\) at each observationi is expressed as:
\(\hat{\beta}\left(u_{i},v_{i}\right)={[H^{T}W\left(u_{i},v_{i}\right)H]}^{-1}H^{T}W\left(u_{i},v_{i}\right)Z\), (5)
where
\(H=\par \begin{bmatrix}1&{h}_{11}\cdots&{h}_{1L}\\ \vdots&\ddots&\vdots\\ 1&{h}_{n1}\cdots&{h}_{\text{nL}}\\ \end{bmatrix}\), \(Z=\par \begin{bmatrix}{z}_{1}\\ \vdots\\ {z}_{n}\\ \end{bmatrix}\), (6)
\(W\left(u_{i},v_{i}\right)=diag(w_{1}\left(i\right),\ \ldots w_{j}\left(i\right)\ldots{,w}_{n}\left(i\right)\)), and (7)
\(\hat{\beta}\left(u_{i},v_{i}\right)={(\beta_{0}\left(u_{i},v_{i}\right),S_{k1}\left(u_{i},v_{i}\right){,\ldots,S}_{\text{kL}}\left(u_{i},v_{i}\right))}^{T}\). (8)
\(H\) and \(Z\) represent the matrix form of \({h}_{\text{il}}\) and\({z}_{i}\) for the n observations and L layers used in this study. \(W\left(u_{i},v_{i}\right)\) is a spatial weight matrix based on the Euclidean and Gaussian distance decay-based functions in the spatial domain. In this case, the Euclidean distance is expressed in meters. The parameter \(w_{j}\left(i\right)\) between observationsi and neighboring j in the spatial weight matrix is commonly used as a kernel and represents a Gaussian distance decay-based function,\(w_{j}\left(i\right)=exp\left(-\frac{D_{\text{ij}}^{2}}{b^{2}}\right)\), where \(D_{\text{ij}}\) is the Euclidian distance between observationsi and j .\(D_{\text{ij}}=\sqrt{{(u_{i}-u_{j})}^{2}+{(v_{i}-v_{j})}^{2}}\)and \(b\) is the non-negative parameter known as the bandwidth (Brunsdon et al., 1996; Chu et al., 2018). Bandwidth b can be determined by several criteria, such as the cross validation (CV) procedure (Fotheringham et al., 2003). The CV procedure is used to select the optimal parameters b :
CV(b ) =\(\sum_{i}{(\text{Δz}_{i}-\text{Δz}_{\neq i}(b))}^{2}\) . (9)
For comparison, the models for linear and spatial regressions i.e., OLS and SR are applied on all four aquifer layers, including mixed and single layer drawdowns from different aquifers. Spatial estimation of the drawdown is applied using inverse distance weighting (IDW). The models are developed using MATLAB. Model 1 is based on the cumulative drawdown from all aquifer layers. Models 2–5 are based on only drawdowns from Layers 1–4 without the y-intercept, respectively, as shown at a later stage.
To evaluate the model performance, the common measures R2 and RMSEs are used on the basis of observation data and estimated values at these points.
Results and discussion
The subsidence–accumulated drawdown relation was built in linear and spatial regressions; that is, the SR model. First, a relation between total subsidence and multi-layer drawdowns was accurately formulated, and the spatial patterns of subsidence were subsequently estimated.