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.