Huang et al. Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-020-00346-y(2020) 7:51

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Identification of hydrogeological evolution
using hydrogeology-seismology analysis ofgroundwater head fluctuation related to
the 1999 MW = 7.5 Chi-Chi earthquake
Chien-Lin Huang, Nien-Sheng Hsu*, Chun-Hao Yao and Gene J.-Y. You

Abstract
At 1:47 a.m. on September 21, 1999, the Mw 7.5 Chi-Chi earthquake struck Taiwan. The purpose of this study is to(1) apply multiple spatiotemporal-frequency analysis to filter the post-seismic change in groundwater head to
explore the implicit drawdown associated with excess pore-water pressure release and effective stress relief during
post-seismic evolution and (2) establish a stochastic experimental optimization model for identifying thehydrogeological evolution. The approaches used for post-seismic drawdown filtration include multi-rank principal
components decomposition, multi-frequent wavelet transforms decomposition, and multi-level wavelet de-noising.
This study especially evaluates the following advanced post-seismic evolving parameters: (1) harmonic averageleaking/injecting rate, (2) distance between the acting position and monitoring well, (3) storage coefficient undereffective stress relief and formation compression, and (4) transmissivity for excess-pore-water pressure release. This
study applies the integrated methodology on 179 monitoring wells in the Chou-Shui River alluvial fan. Results showthat the overlying principal components PCs and low-level wavelet de-noising can filter additional sources/sinks, inwhich the extracted drawdown from PC1+PC3 was related to the excess pore-water pressure relaxation process,that from PC2+PC5+PC6+PC7 and high-frequency wavelet de-noised detail cD2 related to the earth tidal fluctuation
effect, and that from PC4+PC9 and cD3 related to the barometric effect. According to the Riemann integral and an
objective function value duration curve, calculated occurrence probability from the stochastic optimization for SC2,
the storage coefficient was reduced from pre-seismic pumping test value 0.00107, post-seismic 27th hour evolving
value 0.000826 to post-seismic pumping test value 0.000578 in 2004, and the transmissivity increased from preseismic test value 92.4 m2/h, post-seismic 27th hour evolving value 98.6 m2/h to post-seismic test value 147.6 m2/h.
The results demonstrate that the SC2 and GH3 zones suffer from crustal compression and the permeability wasincreased to dissipate excess pore-water pressure and effective stress.
Keywords: Hydrogeological evolution, Hydrogeology-seismology analysis, Co-post-seismic groundwater headfluctuation, Earthquake, Excess pore-water pressure, Crustal compaction* Correspondence: nsshue@ntu.edu.tw
Department of Civil Engineering, National Taiwan University, No. 1, Sec. 4,
Roosevelt Road, Taipei 10617, Taiwan© The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in this article are included in the article's Creative Commonslicence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commonslicence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain
permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

(2020) 7:51

Page 2 of 28

Introduction
The crustal stress (strain) caused by earthquakes is anatural driving force, causing the groundwater head tochange in a short period in a large area subject to a consistent stressor (Wakita 1975). Hence, it provides a
medium for understanding the hydrogeological evolutionof the formation through earthquake-induced changes ingroundwater head. The changes in groundwater headcaused by earthquakes can be divided into co-seismic oscillation and step-like changes (Wang et al. 2004a). The
co-seismic oscillation is mainly subject to the elasticchange caused by seismic wave transmission, in which
the groundwater head returns to the original head after
the seismic wave dissipates. The step-like changes ingroundwater head are mainly affected by the inelasticchange of the crustal stress (or strain), and the groundwater head does not return to the original head immediately after the earthquake, resulting in possible
permanent head changes after a period of time.
The instantaneous rise and fall of groundwater head
induced by the earthquake are mainly produced by thefact that the stratum is squeezed by the seismic stress,
causing the water pressure of the large area to rise, sothat the excess pore-water pressure may be dissipated tothe surface. The co-seismic step-like change in groundwater head is usually accompanied by a post-seismic recovery phenomenon, which represents the pore-waterpressure in the stratum stabilizing. This behavior ismainly controlled by the hydrogeological characteristicsof the porous medium that can be identified from thedissipation or replenishment of pore-water pressure(Roeloffs 1996).
Three categories of previous studies regarding abnormal changes in groundwater head caused by earthquakes
can be classified according to timed changes in groundwater head. In research focusing on pre-seismic hydrogeological abnormalities, Roeloffs and Quilty (1997),
Koizumi et al. (2004), and Hartmann and Levy (2006)
found that some monitoring wells recorded abnormal
groundwater head and hydrochemical changes before
earthquakes. Regarding research focusing on abnormalities in co-seismic groundwater head change, Roeloffs(1996), Wakita (1975), and Kagabu et al. (2020) foundthat the regional distribution of co-seismic groundwaterhead rise and fall closely coincided with areas of contraction and dilatation expected from faulting. In research focusing on post-seismic hydrogeologicalchanges, Rojstaczer and Wolf (1992), Montgomery andManga (2003), and Petitta et al. (2018) found that earthquakes result in persistent or permanent changes inaquifer characteristics, hydrogeological parameters, andgroundwater hydrology.
The 1999 Chi-Chi (Mw 7.5) earthquake was the largest
to hit Taiwan in the past century. A dense network ofhydrological stations on the Chou-Shui River alluvial fan
near the epicenter of the earthquake captured thegroundwater head changes before, during, and after theearthquake. The proximity to a large earthquake and thedense network of hydrological stations provided a rareopportunity to examine the effect of an earthquake ongroundwater. Several studies have discussed the abnormal groundwater head changes in response to this earthquake (Chia et al. 2008; Wang et al. 2016; Wang et al.
2004b). However, large earthquakes may cause partial
aquifer/aquitard rupture and an increase in the permeability of aquifer/aquitard (Elkhoury et al. 2006;
Hasegawa et al. 2019; Katayama et al. 2015; Ujiie andKimura 2014), but few previous studies considered theearthquake-produced leaking/injecting effect on the
post-seismic changes in groundwater head drawdown.
Moreover, the post-seismic pore-elastic mechanical parameters and hydrogeological parameters associated withexcess pore-water pressure release and effective stressrelief should be in evolution, so many parameters may
increase the uncertainties and error of simulating thedrawdown in groundwater head. Furthermore, previous
studies have seldom used a stochastic approach to simulate all possible natural hydrogeological scenarios, inwhich the measured groundwater head may involve additional factors such as sources/sinks and boundary conditions, where the corresponding magnitude may befiltered by proper spatiotemporal-frequency analysis.
Accordingly, the purpose of this study is to (1) apply
multiple spatiotemporal-frequency analysis approaches tofilter the post-seismic change in groundwater head to explore the implicit drawdown associated with excess porewater pressure release and effective stress relief during an
evolving period and (2) establish the stochastic experimental optimization model for identifying the hydrogeological evolution under the leaking/injecting effect. To
identify the hot zones suffering from intense crustal straincaused by an earthquake, this study investigates/detects
the three-dimensional concussive patterns of spatiotemporal change in co-post-seismic groundwater head and obvious permeability increase. The approaches used for theidentification of drawdown related to excess pore-waterpressure release include multi-rank principal components
decomposition, multi-frequent wavelet transform decomposition, and multi-level wavelet de-noising. Moreover,this study establishes an indicator for evaluating the postseismic hydrogeological harmonic average evolving parameters according to the Riemann integral and an objective function value duration curve.

Huang et al. Progress in Earth and Planetary Science

Methodology
Procedures

The research procedures of this study consist offour parts: (1) explicit global hydrogeology-seismology

(2020) 7:51

Page 3 of 28

analysis across multiple zones, aquifers, and aquitards
during excess pore-water pressure release; (2) implicit
regional post-seismic spatiotemporal-frequency analysis
for non-pore-water pressure release-related groundwaterhead filtration; (3) establishment of the stochastic experimental optimization model; and (4) solution of thestochastic optimization model for identifying the evolving process. The flowchart of the methodology is shownin Fig. 1 and the steps are described as follows:
Step 1. According to the measured co-seismic uplift
and post-seismic drawdown in groundwater heads, thisstudy firstly investigates three natural phenomena inhydrogeology-seismology: (1) vertical concussion hotzones and spatial patterns of changes in co-post-seismic
groundwater head, (2) spatiotemporal pattern of postseismic uplift/downward partial derivatives, and (3) detect the hot zones of vertical permeability increase. Accordingly, the three-dimensional hot zones suffering
earthquake-induced intense crustal strain are identified
and selected.
Step 2. Filter the measured groundwater head associated with the seismically induced crustal strain followedby relaxation from excess pore-water pressure. Threedifferent signal decomposition approaches are used to

Huang et al. Progress in Earth and Planetary ScienceFig. 1 Flowchart of the research methodology

(2020) 7:51

Page 4 of 28

evaluate the strain relaxation process, high-low-rank singular value decomposition (SVD), high-low-frequency
wavelet transform decomposition, and high-low-levelsoft-thresholding de-noising. These approaches are usedin conjunction with the measured earth-oceanic tidal
head, barometric pressure, and other geophysical factors.
Step 3-1. Embed a post-seismic groundwater headdrawdown simulation model into an established stochastic experimental optimization loop: (1) set an initial solution for condition ID number i, boundary condition b,
and drawdown analysis scheme α; (2) set the parameterized decision variables (i.e., harmonic average leaking/injecting rate Qð F; r; zÞ, the distance between harmonic
acting position and monitoring well rð F; r; zÞ , storage
coefficient under effective stress compression SðG; η; S s ;F; r; z; SÞ, and transmissivity for excess pore-water pressure release T ð F; r; z; TÞ ; and (3) establish the postseismic groundwater head drawdown simulation modelby transforming the force point source analytical solution of pore elasticity theory (Roeloffs 1996) into poreelastic homologous Theis equation to reduce the number of parameters and uncertainties under crustal straininduced leaking/injecting and excess pore-water pressure
release through the evolving parameters.
Step 3-2. Establish the stochastic optimization model,
which is composed of a stochastic objective function
g α;b;i ðQ; r; S; T Þ , the physical constraints of harmonic
variables Qð F; r; zÞ, rð F; r; zÞ, SðG; η; S s ; F; r; z; SÞ, and T
ð F; r; z; T Þ, in which the objective duration curve (ODC)
of exceedance probability vs. objective function value iscomputed from a series of stochastic experimental calculated gα, b, i.
Step 4. Apply the simulated annealing algorithm tosolve the stochastic optimization model, and output all
possible excess pore-water pressure release scenarios.
Evaluate the post-seismic hydrogeological evolving parameters according to Riemann integral and the computed ODC, and finally discuss the pre-seismic andpost-seismic evolution.

∣≥
vs. time should be larger than 5 cm, i.e., j f∂t
5 cm.
(2) Investigate the concussion change in co-postseismic groundwater head across multiple zones κ∂h ðt;κÞand aquifers f ( f∂t ) by using quartile analysis andthe spatial patterns of change in co-seismic ground-

Huang et al. Progress in Earth and Planetary Science

Global spatiotemporal investigation of hot zones patternof change in co-post-seismic groundwater head anddetection of permeability increase

This study first analyzes the explicit information embedded in the step-like sudden rise/drop in co-seismic
groundwater head followed by a post-seismic release/
recovery behavior. The designed research topics
include:(1) Investigate the hot zone of intense changes in copost-seismic groundwater head h across multiple
zones κ and aquifers f using quartile analysis, inwhich the partial derivatives of groundwater head∂h ðt;κÞ∂h ðt;κÞfwater head across multiple aquifers (i.e., ∂t).
0(3) Investigate the spatiotemporal pattern of postseismic uplift/downward partial derivatives ofgroundwater head vs. time during excess porewater pressure release periods t1 and
t2 across35).
multiple zones and aquifers (Εt ½∂h ∂tðt;κÞ f

 t ¼ t1  t2

(4) Detect the hot zone patterns of obvious verticalpermeability increase during excess pore-waterpressure release evolving periods t1 and t2 across
multiple zones and aquitards T (Εt ½∂h ∂tðt;κÞ  ∂h ∂tðt;κÞ 
f

3

f þ1

 t ¼ t1  t2

5), in which the aquitard T is coated by aquifer f and
f + 1. The condition∂h ðt;κÞ ∂h ðt;κÞ 
Εt ½ f∂t  f þ1∂t 

 t ¼ t1  t2

3

5<0

for the zone κ means that the groundwater head in
one of the aquifers among f and f + 1 rises and thatof the adjacent aquifer falls, indicating a leak to another aquifer through aquitard T.
Hydrogeology-seismology spatiotemporal-frequency
analysis of groundwater head fluctuation to filter factors
related to the strain and relaxation process

According to Roeloffs (1996), the phenomenon of thestep-like sudden rise/drop in co-seismic groundwaterhead followed by a release/recovery behavior is mainly
caused by crustal strain followed by relaxation. However,the post-seismic-measured groundwater heads might beaffected by factors (source/sink) such as rainfall Rrain,
artificial pumping Phuman, earth tidal Tearth, oceanic tidal
Tocean, barometric effect Bbaro, and other factors O(which include well loss, wellbore storage effects, andnoise from automatic recording). These existing sources/
sinks, which are formulated as W in Eq. (1), will affectgroundwater flux. The governing equation of the threedimensional groundwater flow can be expressed as(McDonald and Harbaugh 1988):






∂
∂h
∂
∂h
∂
∂h
∂h
K xx
þ
K yy
þ
K zz
þ W ¼ Ss
∂x
∂x
∂y
∂y
∂z
∂z
∂t

ð1Þ
where Kxx, Kyy, and Kzz are the hydraulic conductivities
along the x, y, and z directions, respectively [LT−1]; Ss isthe specific storage of porous media [L−1]; h is thegroundwater head [L]; and W is the volumetric flux per

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 5 of 28

unit volume representing sources (recharge) and/or
sinks (pumpage) of the groundwater system.
To filter the post-seismic drawdowns that are associated with the crustal strain followed by the relaxation
process of excess pore-water pressure, this study applies
multi-rank principal component analysis, multifrequency wavelet transform, and multi-level wavelet denoising to decompose the groundwater head fluctuation
into a series of intrinsic mode functions (IMFs). The detailed calculation methods are described in the following
sections.

To filter out the factors that are not associated with
post-seismic pressure relaxation among the L factors,
this study uses the partial component j of the entry-wisefactor score ct; j; f , factor loading ϕ j;κ f ; f , and sample
standard deviation for scale 1=εκ f ; f to reduce the nonrelated factors. The post-seismic relaxation processrelated groundwater head hpost‐sei‐relax
and the corret;κ f ; f

Principal component analysis (PCA)− hpost‐sei‐relax
can be gensponding drawdown hpost‐sei‐relaxt;κ f ; f0;κ f ; fLonguevergne et al. (2007) firstly applied PCA to decompose the groundwater head hydrograph into a series of
intrinsic mode functions from low-rank to high-rankcomponents that were assigned to different root causes.
This study uses single or high-low rank principal component combinations to identify the factors/causes associated with the post-seismic pressure relaxation process.
Singular value decomposition (SVD) used in PCA (Pearson 1901) can produce a series of non-correlated newvariables (partially measured groundwater head by a geophysical drawdown factor) for the analyzed spatiotemporal variables (originally measured groundwater head),and the deviation matrix X can be expressed by the fac~ and the factor loading F as follows (Longuetor score C
vergne et al. 2007):
X ¼ UΣVT ¼ UΣVT D‐1 DTpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1~ FT D¼ T s − 1U pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ D‐1 VΣ D ¼ C
Ts − 1

ð2Þ

where the row vector of the n × κf-order matrix U is anorthonormal left singular vector {u1, · , f, …, un, · , f}; D¼ diagðε1; f ; …; εκ f ; f Þ; Σ ¼ diagðσ 1; f ; …; σ κ f ; f Þ, and σ 1; f
≥ … ≥ σ κ f ; f ≥ 0 is a singular value. The row vector of theκf × κf-order matrix V is an orthonormal right singular
vector fv1;; f ; …; v κ f ;; f g . The singular value σj, f deter1 2
mines the eigenvalue λ j; f ¼ n‐1
σ j; f . Next, let F ¼ ½ϕ j;κ f ; f be the κf × κf-order factor-loading matrix, which ϕ j;κ f ; f
is the correlation coefficient between the κ thf monitoring
well’s groundwater head and the jth principal component
coefficient.
The mathematical meaning in Eq. (2) can be interpreted as there being L independent factors causing thedrawdown in groundwater head. The factor score (transpose) matrix CT ¼ ½c1;; f ⋯ cn;; f  stores n time serial
observations for L = κf factors in aquifer f, wherect;; f ¼ ðct;1; f ; …; ct; j; f ; …; ct;L; f ÞT . Hence, Eq. (2) can be
rewritten as:

xt;κ f ; f ¼ ht;κ f ; f − h0;κ f ; f ¼

κf
X



ct; j; f  ϕ j;κ f ; f  1=εκ f ; f

ð3Þ

j¼1

erated as follows:
hpost‐sei‐relax− hpost‐sei‐relaxt;κ f ; f0;κ f ; f
¼ θ1  ct;1; f  ϕ 1;κ f ; f  1=εκ f ; f þ θ2  ct;2; f  ϕ 2;κ f ; f  1=εκ f ; fθ j ∈½0; 1
þ… þ θκ f  cg
t;κ f ; f  ϕ κ f ;κ f ; f  1=εκ f ; f s:t:ð4Þwhere θj is the binary variable to determine whether
component j is related to the post-seismic relaxation
process. If component j is related, then θj = 1; otherwise,θj = 0.
Wavelet transform

The application of wavelet transform to groundwater,hydrology, and water resources are still very new.
Smith et al. (1998), Cannas et al. (2006), and Adamowski and Chan (2011) found that wavelet transform is a feasible data pre-processing method to isolate
or filter relevant factors. Meyer (1990) and Mallat (1989)
used multi-scale analysis to decompose complex signals
through a wave filter into approximation space/coefficients (low-frequency major component) and detail space/coefficients (high-frequency detailed component), as
shown in Fig. 2a.
This study applies the scaling filter associated with theapproximation coefficient and reduces partially leveled
detail coefficient to identify the groundwater head drawdown factors which are associated with the post-seismicpressure relaxation process, where uses the fast discrete
Daubechies wavelet transform (Daubechies 1988) to calculate the above filters and coefficients. This gives a representation in basis functions of the corresponding
subspaces using a linear combination of the scaling function φ(x) and the wavelet function ψ(x), as expressed
below:~f ðxÞ ¼

− γ
2ν X
−1

k¼0

ðν − γ Þ

ak

φ

x
2

γ
γ −2 kþ− lγ 2ν X
−1
X

l¼1

k¼0

ðν − lÞ

dk

ψ



x
− 2l k2lð5Þ(2020) 7:51

Huang et al. Progress in Earth and Planetary Science

Page 6 of 28

Fig. 2 a The classical DWT decomposition (Mallat algorithm). b The multi-level DWT decomposition (level-3)ðν − γÞwith the approximation coefficient ak

and detail co-ðν − lÞ
efficient d k
at kth level:

8
>
>
>
>
<

ðν − γ Þ

¼

ak

N
−1
X

:
ðν − γþ1Þη j að2kþ jÞ mod2ν − γþ2

j¼0

N
−1
X
>
>
ðν − lÞ
ðν − lþ1Þ>d¼
ð − 1Þ jþ1 η j að2kþ2N − 1 − jÞ mod2ν − lþ2
>
: k
j¼0

ð6Þ
where γ = 1, 2, …, ν + 1 is the γ-th recursion step
resulting in the linear combination between the multilevel approximation and detail space, and the raw
groundwater head signal x = X(t) underwent three levelsof the wavelet transform, as shown in Fig. 2b. For N ∈ ℕ,
a Daubechies wavelet ψ(x) and scaling function φ(x) ofclass D-2N is defined as follows:
8−1
X
pﬃﬃﬃ 2N>>>ðÞ≔ð − 1Þk η2N − 1 − k φð2x − k Þ
ψ
x
2
>
<
>
>
>
>
:

k¼0

φ ð xÞ ¼−1
X
pﬃﬃﬃ 2N2

ð7Þηk φð2x − k Þk¼0

where N is the number of Daubechies coefficients ηk,and η0, …, η2N − 1 ∈ ℝ are the constant filter coefficientssatisfying the conditions:
N
−1
X

N
−1
X
1
η2k ¼ pﬃﬃﬃ ¼
η2kþ1
2
k¼0
k¼0

2NX
− 1þ2lk¼2l

ηk ηk − 2l ¼

8
< Z φ ð xÞ ¼ 0

1 if
0 if

ð8Þ

ℝ

for

x∈ℝn½0; 2N − 1φð2x − k Þφð2x − lÞdx ¼ 0 for k≠lð10ÞThis study uses the partial Daubechies scaling father
wavelet associated with the corresponding approximation coefficients in the (ν − γ)th-leveled component calculated drawdown in the stochastic experiment
optimization to evaluate the possible post-seismic hydrogeological evolution.
Wavelet de-noising

The principle of wavelet de-noising is to retain the wavelet coefficient to define the factors/causes associatedwith the post-seismic relaxation process. Suppose that
there are n noisy samples of a function fde(t) and we observe datax ¼ ðh1;κ f ; f − h0;κ f ; f ; …; ht;κ f ; f − h0;κ f ; f ; …;hn;κ f ; f − h0;κ f ; f Þ . The most general de-noising problem
for the noisy groundwater head signal is expressed asfollows Donoho (1995):
Xðt Þ ¼ f de ðt Þ þ σeðt Þ;

t ¼ 1; 2; …; n

ð11Þ

where fde(t) is a de-noised function sampled at equally
spaced points and {e(t)} are independent identical distributed (i.i.d.) in N(0, 1). The de-noising objective is tosuppress the signal noise X(t) and recover fde(t). Our
de
de
de
de
goal is to estimate ^f ¼ ½ ^f ð1Þ; …; ^f ðtÞ; …; ^f ðnÞwith small mean square error, as expressed below:
Minimize^ ðν − l Þ
d
k

2
n 
 de
 ðν − lÞ1X^f de t; d^Ω ^f ; f de ¼
Þ − f de ðt ÞÞ ;kn t¼1

k ¼ 0; 1; …; 2ν − l − 1

ð12Þ
l¼0
;
l≠0

^ ðν − lÞ is the de-noised detail coefficient in level ν
l ¼ 0; 1; …; N − 1where d k− l. The noise is concentrated in the high-frequency
ð9Þ IMFs, so those are removed by using the fast discrete
Daubechies wavelet transform with a soft threshold. The

(2020) 7:51

Huang et al. Progress in Earth and Planetary Science

soft-thresholding rule is usually referred to as wavelet
shrinkage since it reduces the high amplitude detail coef-

Page 7 of 28G∇2 uz þ Gð2η − 1Þðν − lÞ

ficients d k
so that their absolute value is below a certain threshold level, denoted by λ(l, k), which cangenerally be a function of the decomposition/resolution
level l and the index k. This study calculates the de^ ðν − lÞ using soft-thresholding
noised wavelet coefficients dkrule δ soft
λ , which the associated threshold λ(l, k) is calculated according to Donoho (1995), as expressed below:^ ðν − lÞ ¼ δ soft d ðν − lÞ ¼
d
k
λ
k

8
>
>
<

0;



 ðν − l Þ 
d k
 ≤λ

if

ðν − l Þ− λ; if d ðkν − lÞ > λdk>>: ðν − l Þdkþ λ; if d ðkν − lÞ < − λð13Þλðl; k Þ ¼ σ ðl; k Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2 ln ðnÞ

ð14Þ

where σ(l, k) is the standard deviation of the noise signal in level l. The last step is to compute wavelet reconde
structed de-noised signal ^f ðtÞ based on the originalðν − γÞof level ν − γ and theapproximation coefficients ak
ðν − lÞ
^of levels 1 to γ.
modified detail coefficients dk∂ε G 2 ∂uθ ur∂Pþþ−∂r r r ∂θ∂rrFrδ ðr Þδ ðθÞδ ðzÞH ðt Þ ¼ 0r

ð15Þ
G∇2 uθ þ Gð2η − 1Þ
þ



−

K 2
∂ε
∂P Q
∇ P þ þ np β f− δ ðr Þδ ðθÞδ ðzÞδ ðt Þ ¼ 0
γw
∂t
∂t
r

ð18Þ

where r is the horizontal distance from the origin, z is
the depth of the force point source, ε is the volume
∂uz
ur
1 ∂uθ
r
strain (ε ¼ ∂u∂r þ r þ r ∂θ þ ∂z ), G is the shear modulus,η = (1 ‐ ν)/(1 ‐ 2ν), ν is the Poisson’s ratio, γw is the specific weight of water, np is porosity, βf is water compressibility, δ is the Dirac delta function, and H is the
Heaviside step function. If the pumping/injecting water
effect (Q = 0) and the tangential force Fθ are ignored, assuming the forces along horizontal and vertical directions are uniformly applied (Fr = FZ = F), letting Dt ' =R2/4Dt and C0 = 1/8πGηSs, where Ss = K/D and D = T/S
is the hydraulic diffusivity which can be estimated fromthe transmissivity (T) and storage coefficient (S). Aftersimplification, the point source force (F) can determine
the change in groundwater head (Δht) at different
positions:P
C 0 F ðr þ z Þ
Δht ¼
¼
γw
R3

s0 ðt Þ ¼

According to Roeloffs (1996), considering a saturated porous medium in an infinite domain, assuming that themedium is homogeneous isotropic soil and is subjected tothree concentrated forces produced by the earthquakeinduced crustal strain along the coordinate axis (Fr, Fθ,and Fz) and the Q-volume pore-water is injected at t = 0.
When t > 0, the concentrated force continues to act on
the medium. If the media displacement ur, uθ, and uz andthe excess pore-water pressure P are the variables and thecoordinate origin is placed on the force point source, thegoverning equation of the pore elasticity theory can be
expressed in cylindrical coordinates:

þ

ð17Þ

"
erf



1

D0 t 2#
 0 12
4D t
0expð − D t Þ
−
π

ð19Þ

Formulations of co-post-seismic groundwater head
recovery equation
Force point source analytical solution

G∇2 ur þ Gð2η − 1Þ∂ε ∂P F z−þ δ ðr Þδ ðθÞδ ðzÞH ðt Þ ¼ 0
∂z ∂z
r



1 ∂ε G 2 ∂ur uθ
1 ∂P
−
−
þ
r
r ∂θ r r ∂θ
r ∂θFθδ ðr Þδ ðθÞδ ðzÞH ðt Þ ¼ 0
r

ð16Þ

t
X

Δhts

ð20Þ

t s ¼1

where R is the distance between the force point source
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
and monitoring point ( R ¼ r 2 þ z2 ), erf(·) is the error
function, and s ' (t) is the drawdown of the monitoring
wells. The parentheses in Eq. (19) represent the excess
pore-water pressure corresponding to the change ingroundwater head was gradually removed as timeincreases.
However, during co-post-seismic pressure release, thehydrogeological properties in the porous media are
evolving, so accurate estimations of G, η, Ss, F, T, S, andcoefficients of the error function in Eq. (19) for different
periods are difficult and should be time-variant, andshould not equal the pre-seismic values.
Reduction of parameter and uncertainty by usingtransformed pore elasticity-based equation

In reality, when an earthquake or swarm earthquake occurs, the forces along the horizontal and vertical directions are usually nonuniform (Fr ≠ FZ ≠ F) and thetangential force Fθ cannot be ignored, because thesource, depth, and scale of each earthquake are different.
Besides, the seismic point source forces often cause formation fracture, which produces the post-seismic

(2020) 7:51

Page 8 of 28

leaking/injecting Q ≠ 0. Let the harmonic average leaking/injecting rate Q act on a position with a harmonic
average distance of r between monitoring wells. To reduce the number of variables/parameters in Eq. (19)along with uncertainties, this study uses the transformed
Theis equation to simulate the post-seismic groundwaterhead drawdown that finishes with a relaxation process.
For the post-seismic groundwater flow, the approximate
analytical solution of the transformed Theis Equation(Theis 1935) can be expressed as follows:

uncertainty/reliability (occurrence probability) of eachgenerated solution to be quantified by the exceedance
probability of the objective function value duration curve(ODC). The lower objective function value (lower uncertainty) in the ODC would result in a higher exceedance
probability. The filtered measured drawdown that ismore like the actual relaxation process and better fit between the filtered drawdown and calculated drawdown
would result in a lower objective function value. The optimized solutions corresponding to the objective function value is within a set threshold gthresh surrounded
minimum that is sent to evaluate the post-seismic release scenario through hydrogeology. This study uses theroot mean square error (RMSE) as an objective function
gα, b, i to calculate the fitting accuracy between the measured/filtered drawdown s0 mea
α ðtÞ and the calculated
0 cal
drawdown s b;i ðt; Q; r; S; TÞ, as expressed below:

Huang et al. Progress in Earth and Planetary Science∂2 h
1∂h S ðG; η; S s ; F; r; z; S Þ ∂h
2 þ r ð F; r; zÞ ∂r ¼
∂t
∂r
T ð F; r; z; T ÞQð F; r; zÞs ðt Þ ¼ H ‐hðt Þ ¼
4πT ð F; r; z; T ÞQð F; r; zÞ¼
W ðuÞ
4πT ð F; r; z; T Þ
0

u¼

0

Z ∞
0

ð21Þ

e−u
du
u
ð22Þ

Minimizer 2 ð F; r; zÞS
4T ð F; r; z; T Þt

W ð uÞ ¼

Z ∞
0

e

−u

u

du ¼ − 0:5772 − ln u þ u −

ð23Þ
2

3

Q;r;S;T

ð25Þ

4

u
u
u
þ
−
þ⋯
2  2! 3  3! 4  4!

ð24Þ
0

where H is the initial groundwater head, Qð F; r; zÞ isthe harmonic average leaking/injecting rate during a
post-seismic pore-elastic release/recovery period, T ð F; r;z; T Þ is the post-seismic transmissivity of the aquifer forexcess pore-water pressure release with the leaking/injecting effect, rð F; r; zÞ is the harmonic average distance between the leaking/injecting position and monitoring well, SðG; η; S s ; F; r; z; SÞ is the post-seismic
aquifer storage coefficient with effective stress relief andformation compression, and W ðuÞ is the harmonic average well function.
From Roeloffs (1996) to this section, although the transformed Theis equation cannot estimate the groundwaterhead at t = 0, the numbers of variables/parameters for estimating post-seismic recovery groundwater head/drawdown
were reduced from over 8 to 4. However, the uncertainty
still exists that every possible solution along with reliableoccurrence probability should be thoroughly considered.
Establishment of the stochastic experimental optimization
model

This study thoroughly considers all of the possible excess pore-water pressure release-related parameter (Q; r;
S; T ) solutions in a specific post-seismic period under
different kinds of initial conditions i, and boundary conditions b that best fit the filtered drawdown associatedwith the relaxation process using ensemble analysis
scheme ID number α. This study designs thevﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ uXn 
2 α ¼ 1  Α
u
mea
cal
u
s0 α ðt Þ‐s0 b;i ðt; Q; r; S; T Þ  b ¼ 1  Β t t¼1
i ¼ 1  Ι
g α;b;i Q; r; S; T ¼
n

s:t:





g α;b;i Q; r; S; T ≤g thresh ¼ ODC ρthresh
;;

where s0 mea
α ðtÞ is the drawdown associated with the pressure relaxation filtered using analysis scheme ID number
α; s0 calb;i ðt; Q; r; S; T Þ is the calculated drawdown from Eq.
(22) under boundary conditions b and initial condition i;Α, Β, and Ι are the total number of signal analysis
schemes, boundary conditions, and initial conditions, reis the exceedance probability
spectively; and ρthresh
;;
threshold corresponding to the objective function value
threshold gthresh, in which the transformation between
and gthresh is made by the ODC.
ρthresh;;
From Eqs. (15) to (18), one can find that F is positively
proportional to Qr ( F∝ Qr ). Hence, the upper and lower
bounds of variables Qð F; r; zÞ and rð F; r; zÞ of an analyzed monitoring well-controlled zone ν can be evaluated by those of a reference well ρ with post-seismic
pumping test according to the distance r∗ between thecrustal strain-induced fault (leaking/injecting source)
and monitoring well, as expressed below:
8 LB
UB
>
Qρ ð F; r; zÞ ≤Qρ ð F; r; zÞ ≤ Qρ ð F; r; zÞ
>
>
>
LB
UB
>
>
>
< Qρ ð F; r; zÞ ≤Q ð F; r; zÞ ≤ Qρ ð F; r; zÞv
r ν =r ρr ν =r ρ>> r LB ð F; r; zÞ ≤r ð F; r; zÞ ≤r UB ð F; r; zÞ
>
>
ρ
ρ
>
ρ
>
>
: r LB ð F; r; zÞ  r  =r  ≤r ð F; r; zÞ ≤ r UB ð F; r; zÞ;κ∈½ρ;ν
v
ρ
ν ρ
κ

ð26aÞ
ð26bÞ
ð26cÞ
ð26dÞ

where Qρ ð F; r; zÞ and Qν ð F; r; zÞ are the harmonic average leaking/injecting rates of the referenced and nonLB

referenced monitoring zone, respectively; Qρ ð F; r; zÞ(2020) 7:51

Page 9 of 28

and Qρ ð F; r; zÞ are the lower and upper bounds of theharmonic leaking/injecting rate of the referenced monitoring zone, respectively; r ρ and r ν are the distances between the strain-induced fault and monitoring well of
the referenced and analyzed monitoring zone, respectively; r ρ ð F; r; zÞ and r ν ð F; r; zÞ are the harmonic distances between the leaking/injecting position andmonitoring well of the referenced and analyzed monitorUB
ing zone, respectively; and r LB
ρ ð F; r; zÞ and r κ ð F; r; zÞare the lower and upper bounds of the harmonic distance between the leaking/injecting position and monitoring well of the referenced zone ρ and all controlled
zone κ, respectively.
To guarantee the reality, reliability, and accuracy ofthe optimized solution from a series of stochastic experimental optimization, this study sets the upper and lower
bounds for decision variables according to the in situ
soil formation and material, as expressed in Eq. (29). Accordingly, the governing equation of the groundwaterhead drawdown calculation, expressed in Eqs. (22)–(24),should also be included as constraints.
Moreover, in this study, based on the Riemann integral
and mean value theorem, the optimized stochastic situational solutions of each exceedance probability andprobability density were used to calculate the harmonic
hydrogeological evolving parameters S κ ðG; η; S s ; F; r; z; SÞ and T κ ð F; r; z; TÞ . First, the expected post-seismic
evolving parameter of each exceedance probability (EP)in the ODC is imported for the derivation to consider all
scenarios. When calculating the expected parameters ofof the zone κ, respectively; Pκ, ρ is the ρth EP of the evolution parameter; and ρmax is the maximum serial number of EP. The above-formulated optimization problem
is non-linear. This study solved it by using the simulated
annealing algorithm in the MATLAB optimization
toolbox.

Huang et al. Progress in Earth and Planetary Science

UB

E

two adjacent EPs such as S κ ðP κ;ρ  Pκ;ρþ1 ; G; η; S s ; F; r; z
E; SÞ and T κ ðP κ;ρ  Pκ;ρþ1 ; F; r; z; TÞ , the occurrence
probability is the difference between the two EPs Pκ, ρ + 1− Pκ, ρ, while the evolving parameter value is their averS ðP ;ÞþS ðP ;Þ
T ðP ;ÞþT ðP ;Þage, namely ð κ κ;ρ 2 κ κ;ρþ1 Þ and ð κ κ;ρ 2 κ κ;ρþ1 Þ.
Then, to evaluate the mixed influence of pore-waterpressure release, effective stress relief, formation compression, and leaking/injecting effect, the harmonic
evolving parameters are calculated from the summationof the expected parameters of all different EPs, as derived below:
8ρmax  >XSκ P κ;ρ ;  þ Sκ P κ;ρþ1 ; 
>
HAE
>
>
 P κ;ρþ1 − P κ;ρ
S
ð

Þ
¼
>
< κ
2
ρ¼1




ρmax >X
T κ P κ;ρ ;  þ T κ P κ;ρþ1 ; 
>
HAE
>
>
 P κ;ρþ1 − P κ;ρ
T
ð

Þ
¼
> κ
:
2
ρ¼1

ð27Þ
HAE

HAE

where S κ ðÞ and T κ ðÞ are the post-seismic harmonic average storage coefficients under effective stresscompression and the transmissivity for pressure release

Solving stochastic optimization model using simulated
annealing algorithm

This study used a simulated annealing algorithm proposed by Kirkpatrick et al. (1983) to optimize/solve theminimum value of the objective function for different
initial solutions i, boundary conditions b, and analyticalschemes α. The decision variables are four excess porewater pressure release-related parameter solutions considering the leaking/injecting rate Qð F; r; zÞ, acting distance rð F; r; zÞ between the position of Q and amonitoring well, storage coefficient under effective stressrelief-compression SðG; η; S s ; F; r; z; SÞ, and transmissivity for excess pore-water pressure release Tð F; r; z; T Þ .
The optimization flowchart using simulated annealing isshown in Fig. 3, and the optimization steps are described
below:
Step 1. Set initial solution i, boundary condition b, andanalysis scheme α, and then set the initial temperature0

T 0c ¼ 100 and initial decision variables ξ0 ¼ ½Q; r; S; T  .
Accordingly, input the settings into the post-seismic
pore-elastic groundwater head drawdown simulation
model to acquire the analytical solution and calculate
the initial objective function value g0.
Step 2-1. Randomly generate a neighborhood decision0 ττvariable ðξ Þ ¼ ½Q0 ; r 0 ; S 0 ; T 0  and input it into the postseismic groundwater head drawdown simulation model
to calculate the analytical solution, then calculate thecorresponding objective function value (g')τ for a searching iteration ID number τ.
Step 2-2. If (g')τ ≤ gτ, then accept the neighborhood solution, ξτ + 1 = (ξ')τ. If (g')τ > gτ, then a number ϑτ that lies
between 0 and 1, 0 < ϑτ < 1, is randomly generated andthe Metropolis acceptance criterion (Metropolis et al.
1953) is used to calculate the acceptance probability, i.e.,0 ττP τ ¼ 1= expððg ÞT τ− g Þ , wherecT τc

is the current

temperature at iteration τ associated with the coolingrate (i.e., T τc ¼ T 0c  0:95τ ). The annealing parameter τ isthe same as the iteration number until re-annealing. Whenϑτ < Pτ, the new variable (ξ')τ is accepted, ξτ + 1 = (ξ')τ; otherwise, it is not accepted, and the process returns to step 2-1.
Step 3. Examine whether ξτ + 1 meets the stopping condition (i.e., the gradient of objective function gτ vs. iter∂g
is lower than 1 × 10−10 after 2000
ation times ∂τiterations). If the stopping condition is satisfied, thenoutput the optimized g value and the correspondingHuang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 10 of 28

Fig. 3 The flowchart of identifying post-seismic hydrogeological evolution using simulated annealing algorithm

decision variables. If the condition is not satisfied, a reannealing scheme is added to strengthen the ability for
identifying the global minimum. If 100 more new variables during the past iterations are accepted, the annealing
parameterτisreplaced
by

1


∂g
A
max ∂τ 

τ
τ
¼
1

100
T

τ ¼ τ 0 ¼ log T 0τ
to increase

c

∂g 
∂τ 
 τ ¼ 100
τ0
T c for jumping out of local minima; otherwise, if it has
not accepted, then set τ = τ + 1 and return to step 2-1.

!

Application
Study area and pre-seismic hydrogeology

The study area of this research was the groundwater system of the Chou-Shui River alluvial fan, located in
middle-west Taiwan. This alluvial fan is mainly

recharged by the Chou-Shui River, which extends fromthe Baguah Mountain and Dulliu Hill in the east to the
Taiwan Straits in the west and from Wu River in thenorth to Peikang River in the south. The radius of the alluvial fan is approximately 40 km and the area is approximately 2079 km2. The Chou-Shui River alluvial fan
can be divided into proximal, mid, and distal fan regions
according to topography, geology, and strata, as shownin Fig. 4a.
The groundwater strata in the Chou-Shui River alluvial fan is approximately 330 m thick, which can bedivided into four aquifers and four aquitards, as
shown in Fig. 4b. Aquifer 1 covers the whole region
and varies in thickness from 19 to 103 m. The proximal fan is mainly composed of gravel and coarse
sand, which grades into fine sand and clay in the
middle fan and distal fan. The depth of aquifer 2 is
35 to 217 m from the ground surface, and pre-seismic
S is about 1.15 × 10−4 to 2.98 × 10−3 (Jang et al.
2008). Aquifer 3 covers the region at a depth of 140

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 11 of 28

Fig. 4 (a) Study area. (b) Schematic diagram of Chou-Shui River alluvial fan groundwater strata and evolution comparison of pre-seismic andpost-seismic groundwater head in aquifers 1 (f-1), 2 (f-2), and 3 (f-3)

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 12 of 28

Fig. 5 a Measured co-post-seismic groundwater head fluctuation in Chou-Shui River alluvial fan during the 1999 MW = 7.5 Chi-Chi Earthquake. bComparison of pre-co-post-seismic change in groundwater head

to 275 m, and S is about 5.66 × 10−5 to 1.07 × 10−3.
The pre-seismic average transmissivity (T) of confined
aquifers 2 and 3 is 50.63 and 29.59 m2/h, respectively.
The comparison of pre-seismic and post-seismic
groundwater head along profile A-B is shown in Fig.
4f, in which only well Tanaka-2(TN2) showed dropping in co-seismic groundwater head and the others
showed rising.

Measured hydrogeology-seismology nature resulting
from the strain and relaxation process

At 1:47 a.m. on September 21, 1999, the Mw 7.5 ChiChi earthquake struck Taiwan. The epicenter of thisearthquake was located at 23.85 °N, 120.82 °E, whichis in Chi-Chi Town in Nantou County, at a depth ofapproximately 8.0 km. 157 monitoring wells on thewest side of the Changhua Fault showed oscillatory

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 13 of 28

Fig. 6 Hot zones of intense change in co-post-seismic hourly groundwater head by using quartile analysis in a proximal fan and middle fan andb distal fan

and step-like co-seismic changes in groundwater head,followed by stable drawdown-like decreases/increases
in heads after suddenly rising/dropping, as shown in
Fig. 5a.
Jang et al. (2008) had previously conducted on-site
pumping tests for the SC2, Shihu-2(SH2), and Yanlin2(YL2) monitoring wells in 2004, where only SC2 inaquifer 2 experience the aforementioned post-seismic response. To evaluate the hydrogeological scenario in different aquifers, the GH3 monitoring well in aquifer 3

was also selected as a study subject. The monitoring
wells in the middle fan that experience the seismic response were Yanlin-1 (YL1), Huatan-2 (HT2), andGanghou-1 (GH1) wells in aquifer 2 and Shihu-3 (SH3),
Hesing-2 (HS2), Jiulong-3 (JL3), and Dongguang-4(DG4) in aquifer 3. These seven monitoring wells werecombined with the SC2 and GH3 for a total of ninewells, which the datasets for comparison and subsequent
signal analyses of the pre-co-post-seismic change ingroundwater head is illustrated in Fig. 5b.

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 14 of 28

Fig. 7 Quartile patterns of change in co-post-seismic groundwater head (Δt = 1 h) across multiple sub-fans in a proximal and middle fans and b
distal fan

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Fig. 8 Spatial changing patterns in co-seismic groundwater head across a aquifer 1, b aquifer 2, c aquifer 3, and d aquifer 4

Page 15 of 28

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 16 of 28

Fig. 9 The investigated pattern of post-seismic uplift/downward partial derivatives of groundwater head vs. time during t = 1 and t = 27 acrossmultiple aquifers

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 17 of 28

Fig. 10 The investigated pattern of post-seismic uplift/downward partial derivatives of groundwater head vs. time during t = 28 and t = 55
across multiple aquifers

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 18 of 28

Fig. 11 Detected hot zones of obvious vertical permeability increase across multiple aquitards during a t = 1 and t = 27 h and b t = 28 and t =55 across (-1) aquitard 1, (-2) aquitard 2, and (-3) aquitard 3

Results and discussions
Investigated hot zones and pattern of change in co-postseismic groundwater head across multiple sub-fans andaquifers

Initially, this study investigated the hot zones of intensechange in co-post-seismic hourly groundwater head acrossmultiple sub-fans through quartile analysis, in which thepartial derivatives of groundwater head vs. time (Δt = 1 h)∂h ðt;κÞ

followed j f∂t j ≥ 5 cm , as shown in Fig. 6. From Fig.
6a, this study discovers that the concussive ranges of postseismic groundwater head of the southern upstream anddownstream middle fan are larger than those of the northern middle fan, respectively. This is because the northern
upstream boundary comprises the Changhua Fault, andthe Chiuchiungkeng Fault is only minorly connected withthe southern upstream boundary. Hence, the earthquakeinduced strain force from the Chi-Chi earthquake epicen-

ter was not completely blocked by a fault in the southern
alluvial fan.
The proximal fan suffers huge crustal strain force,so some extreme values of co-seismic groundwaterhead suddenly rise over 60 m. Moreover, Fig. 6b exhibits plenty of extreme co-seismic groundwater heads
in the middle distal fan, indicating that this area experiences large stratum-heightening-induced excess
pore-water pressure. Figure 5b verifies that thestratum was heightened 0.81 to 4.52 m calculated
based on the pre-seismic head (t = 0 h) and postseismic approaching steady head (t = 55 h). The
middle-north distal fan only sees an instantaneous
drop in head, meaning that this area is in the crustal
extension zone, and the middle-south distal fan andproximal fan shows both increase and decrease in coseismic groundwater head.

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 19 of 28

Fig. 12 The decomposed post-seismic drawdown hydrograph in groundwater head across multiple principal components in a SC2 and b GH3
monitoring zone by using SVD

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 20 of 28

Fig. 13 Comparison of filtered drawdown hydrographs of post-seismic groundwater head by using a SVD, b wavelet transform, and c wavelet
de-noising for (-1) SC2 and (-2) GH3

This study investigated the quartile change patterns∂h ðt;κÞ( f∂t )

in co-post-seismic hourly groundwater headacross multiple sub-fans. Figure 7 indicates that theextreme values occur substantially across the proximal, middle, and distal fans, in which the rapid riseand drop in co-seismic groundwater head range from
0.34 to 6.55 m and − 0.31 to − 2.47 m, respectively.
Only the most southern distal fan did not experience
the drop in head. Next, we investigated the spatial
patterns in the co-seismic head changes. Figure 8

shows that the increase is mainly concentrated in the
middle fan, and the decrease dominates at the southern boundary and the northern proximal fan acrossmultiple aquifers.
Investigated spatiotemporal pattern of post-seismic
uplift/downward partial derivatives of groundwater headvs. time

This study investigated the spatiotemporal pattern of
post-seismic uplift/downward partial derivatives ofHuang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 21 of 28

Fig. 14 Comparison of the post-seismic measured and decomposed groundwater head recessional hydrographs by using the a wavelet
transform and b wavelet de-noising for (-1) SC2 and (-2) GH3

3∂h ðt;κÞ5) durgroundwater head vs. time (Εt ½ f∂t 
 t ¼ t1  t2
ing t = 1 and t = 27 across multiple aquifers f. From Fig.
9, this study discovers that in shallow aquifers 1 and 2,
the excess pore-water pressure release rate is faster thanthat in deep aquifers 3 and 4, in which the middle fan in
aquifer 1 and the proximal fan in aquifer 2 already exhibit drawdown-like drops in groundwater head during t= 1 and t = 27. Otherwise, the crustal strain-induced
excess-pore-water pressure is continuously pressurized
at aquifers 3 and 4, except for the southern boundary.

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 22 of 28

Fig. 15 Identified hydrogeological evolution pattern 27 h post-earthquake in stochastic experimental optimization for zones a Si-Chou(2) and b
Gang-Hou(3) associated with the (c) objective function value duration curves (ODCs)Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 23 of 28

Fig. 16 Comparison of the measured, filtered, and simulated post-seismic drawdown hydrograph: (a) SC2 and (b) GH3

Additionally, from Fig. 10, this study demonstrates
that the partial derivative groundwater head in the continuously pressurized zones in Fig. 9 during post-seismict = 28 and t = 55 are transformed into decreasing zones
starting to release the excess pressure. Only some portions of the southern boundary in aquifers 1–3 andnorthern boundary in aquifer 4 continuously suffer
crustal strain pressure.

Detected spatiotemporal hot zones of obvious verticalpermeability increase

Hot zones of obvious vertical permeability increase
acrossmultiple aquitards
and
zones
κ
(
3

∂h ðt;κÞ ∂h ðt;κÞ 5 ) are investigated byΕt ½ f∂t  f þ1∂t 
 t ¼ t1  t2
using the Kriging gridding method. Figure 11 shows that3∂h ðt;κÞ ∂h ðt;κÞ
5 < 0 isthe condition Εt ½ f∂t  f þ1∂t 
 t ¼ t1  t2
met in part of the middle and proximal fans in aquitards 1
and 2 and at the distal fan in aquitard 3 during post-seismict = 1 and t = 27. During the evolving period t = 28 to t = 55,
the temporal rise/drop in groundwater head between the adjacent aquifers which cladding aquitards 1 and 3 in the middle fan and proximal fan all turn into the opposite, so itmeans that the permeability became increasing after t = 28.
Although aquitard 3 did not show permeability increase especially in the proximal fan during t = 1 and t = 27 hours,
the permeability became increasing in proximal fan, middle
fan, and most distal fan during t = 28 and t = 55.

Decomposed pattern of measured co-post-seismic
groundwater heads fluctuation using PCA

From the above four explicit analysis for the step-likesudden rise/drop in co-seismic groundwater headfollowed by a post-seismic release/recovery behavior, it
is clear that the middle fan of the Chou-Shui River alluvial fan is the complex critical evolving pore-water pressure area, which is consistent with the selected
groundwater monitoring wells for hydrogeological implicit evolving scenarios analysis (YL1, HT2, GH1, SH3,
HS2, JL3, DG4, SC2, and GH3). Furthermore, there was
no rainfall, Rrain, at the Chou-Shui River alluvial fan before or after the Chi-Chi earthquake, and the earthquake
resulted in a large blackout, which eliminated possible
impacts from artificial pumping, Phuman.
SVD in PCA was performed on the post-seismic
groundwater head recovery hydrograph for nine monitoring wells. Accordingly, there were nine principal components (nine factors) produced for the analyzed zones
SC2 and GH3, as shown in Figs. 12 and 13a. Results
found that principal component 1 (PC1) and the first
three PCs could explain 97.41% and 99.36% of the originally measured groundwater head, respectively, according to the eigenvalue of each component. From
Figs. 12 and 13a, we show that the decomposed drawdown in groundwater head in PC1, PC2, and PC3 appears to be the refined drawdown and excess pore-waterpressure release/recovery behavior, and the other highrank decomposed drawdown may be related to othergeophysical sources/sinks (earth tidal, barometric pressure, recharge) and noise. Therefore, this study used fourcombinations of PC1, PC1+PC2, PC1+PC3, and PC1+

(2020) 7:51

Page 24 of 28

PC2+PC3 to decompose and filter the drawdown hydrograph from the originally measured groundwater heads,
respectively, to extract the factors associated with the excess pore-water pressure relaxation process, as shown in
Fig. 13a. The filtered drawdown is set as s0 mea
in Eq. (25)αfor optimization.

pumping test for the hydrogeological parameter in 2004;
thus, SC2 was set as the reference well to estimate
the boundary condition of GH3 which did not have apost-seismic pumping test. From Fig. 4, we see theconnection between the epicenter and SC2 monitoring well and between the epicenter and GH3 well all
pass through the Chelungpu Fault outside of theChou-Shui River alluvial fan and Changhua Fault.
This study uses the Changhua Fault as a reference tocalculate the distance ratio r GH3 =r SC2 ¼ 24947 km=818
2km ¼ 3:049 , and then the boundary condition of
GH3 well under the decreasing crustal strain transmissive force according to the distance between thewell and fault, as expressed below:

Huang et al. Progress in Earth and Planetary Science

Wavelet transformed pattern of co-post-seismic
groundwater heads fluctuation

One to three levels of wavelet transform were executedon the post-seismic-measured groundwater heads of
SC2 and GH3 monitoring wells to calculate the approximation coefficients cA1, cA2, and cA3 (low-frequency).
These coefficients, as shown in Fig. 14a, were used to refine the drawdown identifying the composition factors
associated with the relaxation process. The decomposed
drawdown hydrographs using multi-level wavelet transforms are shown in Fig. 13b, in which the decomposed
drawdown is set as s0 mea
in Eq. (25) for optimization.
α
The abovementioned figures demonstrate that thedecomposed low-level drawdown hydrograph shows
more refined drawdown than the high-level drawdownhydrograph, because the Daubechies wavelet filter causes
the high-level drawdown hydrograph at log(t) > 5 h to
exhibit wave shapes.
Wavelet de-noising pattern of co-post-seismic
groundwater head fluctuation

One to three levels of wavelet de-noising were performedon the post-seismic-measured groundwater heads of SC2
and GH3 monitoring wells to eliminate noise and retain
the factors associated with excess pore-water pressure relaxation, as shown in Fig. 14b, and the corresponding
drawdown was calculated as in Fig. 13c. The abovementioned figures show that the de-noised low-level drawdown hydrograph represents the more refined drawdown
and excess pore-water pressure release/recovery behavior
than the high-level de-noised hydrograph, because theapplied soft-thresholding proposed by Donoho (1995)
allows some detail factors in the high-level de-noised
drawdown to be removed. In particular, the multilevel de-noised drawdown hydrographs of monitoring
well GH3 are almost the same as the originally measured drawdown, suggesting that noise in the deep
aquifer is rare.
Identified hydrogeological evolving status at 27 h of postseismic strain and relaxation
Parameter setting of the stochastic experimental
optimization

Before the stochastic optimization, the boundary condition and initial condition should be set properly. The
monitoring well-controlled zone SC2 has a post-seismic

 3 
 3 
150
 m =hr ≤QSC2 ð F; r; zÞ ≤1700 m =hr ; 90m ≤r SC2 ð F; r; zÞ ≤ 2200m
49:196 m3 =hr ≤QGH3 ð F; r; zÞ≤ 557:558 m3 =hr ; 274:411m ≤r GH3 ð F; r; zÞ ≤ 2200m

ð28Þ
Moreover, according to the in situ hydrogeological soil
composition of SC2 and GH3 in a confined aquifer, thecorresponding boundary condition is set as expressed
below:
0:000012 ≤S SC2 ðG; η; Ss ; F; r; z; SÞ ≤0:005; 20 m2 =hr ≤T SC2 ð F; r; z; T Þ ≤ 170 m2 =hr
0:00001≤ S GH3 ðG; η; Ss ; F; r; z; SÞ ≤0:003; 0 m2 =hr ≤T GH3 ð F; r; z; T Þ ≤200 m2 =hr

ð29Þ
Furthermore, this study designs a series of experimentsfor optimization to ensure identifying the global minimum along with an occurrence probability calculated by
the objective duration curve for the objective function
value gα, b, i from different kinds of initial conditions i,boundary conditions b, and analysis schemes α in Eq.
(25). In the solution space (boundary conditions) for theτdecision variables ξτ ¼ ½Q; r; S; T  , Eqs. (28) and (29),this study uniformly generated 4800 combinations of initial decision variables (initial conditions), ξτα;b;i (i = 1,2,…,4800), for the experimental optimization.
Ensemble optimized pattern of pore-water pressure andeffective stress delivered through hydrogeology in SiChou(2)After the stochastic experimental optimization, theevaluated evolving parameters (leaking/injecting rateQð F; r; zÞ, harmonic average distance between the acting position and monitoring well rð F; r; zÞ, storage coefficientundereffective
stresscompression
SðG; η; S s ; F; r; z; SÞ , and transmissivity for excess
pore-water pressure release T ð F; r; z; T Þ) for different
initial conditions i, boundary conditions b, and analysis schemes α were collected. The stochastic optimized solutions forSðG; η; S s ; F; r; z; SÞandT ð F; r; z; T Þ associated with the corresponding objective function value gα, b, i for SC2 is shown in Fig.

Huang et al. Progress in Earth and Planetary Science(2020) 7:51

Page 25 of 28

Fig. 17 The decomposed post-seismic measured groundwater heads for zone SC2 by different filtration analysis methods: a PCA and b level-3
wavelet transform comparing to the measured geophysical data from c ML tide station and d TC barometer station. (Note: The coordinates forthe solid lines are on the left, whereas the coordinates for the dotted lines are on the right.)

15a, in which few local minima of the objective function corresponding to less uncertainly (larger reliability) and higher occurrence probability are revealed.
According to the Riemann integral and an ODC calculated occurrence probability for SC2, as shown in
Fig. 15c, the evolving storage coefficient was reduced
from pre-seismic on-site pumping test value 0.00107,
and post-seismic 27th hour value 0.000826 to postseismic pumping test value 0.000578 in 2004 (Janget al. 2008), and the evolving transmissivity increased
from a pre-seismic value of 92.4 m2/h, post-seismic
27th hour value 98.6 m2/h to post-seismic value of
147.6 m2/h in 2004, indicating that the SC2 zone

experienced large crustal compression, and permeability was increased to dissipate excess pore-water pressure and effective stress. The magnitude of change
between the post-seismic 27th hour evolving value
and post-seismic pumping test value in 2004 showsthat S decreased by 49.6%, whereas T increased by
11.2%. Among all of the stochastic optimized solutions, when the decomposed combination of PC1+
PC3 was used to calculate the drawdown associatedwith the excess pore-water pressure relaxation
process, as shown in Fig. 16a, the optimized objectivefunction value (RMSE 0.06327 m) was the lowest
among the three filtration analysis methods.

(2020) 7:51

Page 26 of 28

Ensemble optimized pattern of pore-water pressure andeffective stress delivered through hydrogeology in GangHou(3)compression zone, its groundwater heads will rise; conversely, if the aquifer is located at a crustal expansionzone, its groundwater heads will fall. Lee et al. (2002)
found that during the Chi-Chi earthquake, the westward
thrusting of the Chelungpu Fault resulted in the westward displacement of the western Changhua Fault.
Therefore, the region west of Changhua Fault is a crustal
compression zone, and thus its groundwater head rises.
Because the SC2 and GH3 monitoring wells bothshowed a rise in co-seismic groundwater head and a fallin the post-seismic identified evolving storage coefficientSðG; η; S s ; F; r; z; SÞ, it is proven that the earthquake resulted in crustal contraction in this region. Moreover,according to Wang et al. (2001), because of the asymmetrical shear stress effects on the soil after an earthquake, when the shear stress is sufficiently large todamage soil, it will expand and crack and resulting in a
post-seismic transmissivity T ð F; r; z; T Þ increase. Meanwhile, Montgomery and Manga (2003) noted that theoccurrence of large earthquakes will cause irreversible
changes in hydrogeological parameters such as storage
coefficient S, hydraulic conductivity K, transmissivity T,
and porosity η. Examination of the results of this study
shows that the trends of the identified harmonic averages for the evolving hydrogeological parameters ofzones SC2 and GH3 conform to the change mechanisms
in post-seismic S and T proposed by the aforementioned
researchers. Therefore, the analyzed post-seismic filtered
drawdown hydrograph associated with the applied stochastic experimental optimization model is a useful
methodology.
By comparing the frequency significance and variation
pattern of the measured geophysical factors (earthoceanic tidal head, barometric pressure) with the decomposed drawdown, this study finds that the extracted
drawdown from PC1+PC3 shown in Fig. 16a could beunderstood as the effects of the excess pore-water pressure relaxation process. According to Bredehoeft (1967)
and Roeloffs (1988), the M2 (12.4 h cycle) tidal constituent was used as the tidal force to calculate its effects on
groundwater heads, and the results showed that a variation of 1–2 cm and 0.4–8 cm in amplitude could be
produced. Similarly, the tide along the seashore near the
Chou-Shui River alluvial fan is also mainly governed by
the M2 (12.4 h cycle) tidal constituent so that there are
two high tides and two low tides each day, in which thetide data from the Mailiao (ML) station measured by theCentral Meteorological Bureau during the Chi-Chi
earthquake is shown in Fig. 17c. The planimetric positions of the ML tide station and TC barometer station
are shown in Fig. 4. Accordingly, the drawdown extracted from PC2+PC5+PC6+PC7 in Fig. 17a and cD2 in
Fig. 17b show the earth tidal fluctuation effect, Tearth.

Huang et al. Progress in Earth and Planetary Science

The stochastic experimental optimized solutions for SðG; η; S s ; F; r; z; SÞ and Tð F; r; z; T Þ associated with thecorresponding objective function value gα, b, i for zone
GH3 is shown in Fig. 15b, in which the entire local minimum shows that the evolving transmissivity increased torelease excess pore-water pressure. According to Riemann integral and an ODC calculated occurrence probability for GH3, as shown in Fig. 15c, the evolving
storage coefficient was reduced from a pre-seismic
pumping test value of 0.000149 to post-seismic 27thhour value 0.000112, which is a decrease of 24.8%; andthe evolving transmissivity increased from a pre-seismic
value of 28.8 m2/h to post-seismic 27th hour value 120.7
m2/h, meaning that the GH3 zone also experienced large
crustal compression and permeability was increased todissipate excess pore-water pressure and effective stressin the deep confined aquifer. Among all of the stochastic
experimental optimized solutions, when post-seismicmeasured groundwater heads were used to calculate
drawdown and optimized parameters, the objective function value (RMSE 0.01221 m) was smaller than that fromother filtration analysis methods. Furthermore, Fig. 14b2 could found that after level 1 to 3 wavelet de-noising
was used to refine the post-seismic measured groundwater head of GH3, the results were identical to those of
the original post-seismic measured groundwater head.
Therefore, it could be deduced that the post-seismic
measured GH3’s drawdown in deep aquifer is only
caused by the excess pore-water-pressure relaxation
process.
Comprehensive discussion

Freeze and Cherry (1979) found that storage coefficient
S is related to the compressibility of water and soil, sothat the S can be expressed as Eq. (30). Accordingly, it
can be seen that when the earthquake compresses theaquifer, the soil porosity η and the coefficient of compressibility of the soil α decrease, so that S also
decreases.
S ¼ ρgbðα þ ηβÞð30Þ

where ρ is the density of water, g is the gravitational acceleration, b is the thickness of the aquifer, and α and β
are the coefficients of soil and water compressibility,
respectively.
According to Roeloffs (1996), the earthquake-induced
instantaneous rise or fall in groundwater heads is usuallyreflected by volume changes caused by the sliding (ordislocation) of aquifers composed of poroelastic media
at seismic faults. If an aquifer is located at a crustal

(2020) 7:51

Page 27 of 28

Moreover, Jacob (1940) found that the changes ingroundwater head and barometric pressure show opposite trends. The barometric data from Taichung (TC) station operated by the Central Meteorological Bureau isshown in Fig. 17d. Accordingly, the drawdown extracted
from PC4+PC9 in Fig. 17a and cD3 in Fig. 17b could beapproximately regarded as the barometric effect, B. With
regard to the drawdown extracted from PC8 in Fig. 17aand cD1 in Fig. 17b, because its magnitude is a cyclic
small hydrograph with periodicity, it could be regarded
as white noise and other interference factors, O.

northern proximal fan in aquifers 1, 2 and 3; (2) compared to the middle fan in aquifer 1 and the proximal
fan, the groundwater head in the continuously pressurized zones at aquifers 3 and 4 except for the southernboundary during post-seismic t = 28 and t = 55 are decreases and starts to release excess pressure with a lagtime of 27 h; and (3) the permeability increased in partof the middle and proximal fans in aquitards 1 and 2
and in the distal fan in aquitard 3 during post-seismic t= 1 and t = 27 h. During t = 28 and t = 55, the middle
and proximal fans in aquitards 1 and 3 all show increasing permeability.
Furthermore, the overlying PCs and low-level wavelet
de-noising can filter the additional source/sink factors,
which affects the fluctuation of the post-seismic drawdown hydrograph. By comparing the frequency significance and variation pattern of the measured geophysicalfactors with the decomposed drawdown, this study
shows that the extracted drawdown from PC1+PC3
could be approximately regarded as the effects of the excess pore-water pressure relaxation process. Accordingly, the drawdown extracted from PC2+PC5+PC6+
PC7 and high-frequency wavelet de-noised detail cD2
could be approximately regarded as the earth tidal fluctuation effect. Moreover, the drawdown from PC4+PC9
and cD3 could be understood as the barometric effect,
and the drawdown from PC8 and cD1 as white noise andother interference factors.
After the stochastic experimental optimization
process, according to the Riemann integral and an
ODC calculated occurrence probability for SC2, thestorage coefficient was reduced from the pre-seismic
on-site pumping test value 0.00107, post-seismic 27thhour evolving value 0.000826, to the post-seismic
pumping test value 0.000578 in 2004, and transmissivity was increased from the pre-seismic value 92.4
m2/h, post-seismic 27th hour value 98.6 m2/h, topost-seismic value 147.6 m2/h. For GH3, the storage
coefficient was reduced from the pre-seismic value
0.000149 to post-seismic 27th hour evolving value
0.000112, which is a decrease of 24.8%, and the transmissivity increased from pre-seismic value of 28.8 m2/
h to post-seismic 27th hour value of 120.7 m2/h.
These results demonstrate that zones SC2 and GH3
suffered large crustal compression resulting in permeability increase to dissipate excess pore-water pressureand effective stress in the deep confined aquifer.

Huang et al. Progress in Earth and Planetary Science

Conclusions
The purpose of this study was to (1) apply multiple
spatiotemporal-frequency analysis approaches to filter
the post-seismic change in groundwater head to explore
the implicit drawdown associated with excess pore-waterpressure release and effective stress relief during postseismic evolution and to: (2) establish the stochastic experimental optimization model for identifying the hydrogeological evolution. To identify the hot zones
experiencing intense crustal strain caused by an earthquake, this study investigates/detects the threedimensional concussive patterns of spatiotemporalchange in co-post-seismic groundwater head and verticalpermeability increase. The signal analysis used to describe pore pressure release during drawdown includes
multi-rank principal components analysis using SVD,
multi-frequent wavelet transform decomposition andmulti-level wavelet de-noising. This study especially
identifies the following advanced post-seismic evolving
parameters: (1) leaking/injecting rate, (2) harmonic average distance between the acting position and monitoring
well, (3) storage coefficient under effective stress relief
and formation compression, and (4) transmissivity forexcess pore-water pressure release related to leaking/injecting effect. A transformed post-seismic drawdown
simulation model based on pore elasticity theory for parameter and uncertainty reduction is embedded into thestochastic optimization model. Moreover, this study establishes an indicator for evaluating the post-seismichydrogeological harmonic average evolving parameters
according to the Riemann integral and an objective duration curve.
At 1:47 a.m. on September 21, 1999, the Mw 7.5 ChiChi earthquake struck Taiwan. This study applied theintegrated methodology to 179 monitoring wells in the
Chou-Shui River alluvial fan. The explicit hydrogeologyseismology analysis results show that (1) the rapid increase in co-post-seismic groundwater head ranged from
0.34 m to 6.55 m and the rapid decrease from − 0.31 to− 2.47 m; where the rising change in co-seismic groundwater head was mainly concentrated in the middle fan,and the decreases were at the south boundary and

Abbreviations
CWT: Continuous wavelet transform; D: Hydraulic diffusivity; DWT: Discrete
wavelet transform; EP: Exceedance probability; F: Point source force;
IMF: Intrinsic mode function; ODC: Objective function value duration curve;
PC: Principal component; PCA: Principal component analysis; RMSE: Rootmean square error; S: Storage coefficient; SVD: Singular value decomposition;
T: Transmissivity

(2020) 7:51

Page 28 of 28

Acknowledgments
We thank Mary who works for the company Editage for the English
language editing services, and we thank the Water Resources Agency,
MOEA, for providing valuable hydrogeology-seismology data. Besides, the authors are indebted to the reviewers for their valuable comments andsuggestions.

Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optimization by simulated annealing.
Science 220:671–680
Koizumi N, Kitagawa Y, Matsumoto N, Takahashi M, Sato T, Kamigaichi O,
Nakamura K (2004) Preseismic groundwater level changes induced by crustal
deformations related to earthquake swarms off the east coast of Izu
Peninsula, Japan. Geophys Res Lett 31
Lee M, Liu TK, Ma KF, Chang YM (2002) Coseismic hydrological changes
associated with dislocation of the September 21, 1999 Chichi earthquake,
Taiwan. Geophys Res Lett 29
Longuevergne L, Florsch N, Elsass P (2007) Extracting coherent regional
information from local measurements with Karhunen-Loeve transform: casestudy of an alluvial aquifer (Rhine valley, France and Germany). Water Resour
Res 43
Mallat SG (1989) A theory for multiresolution signal decomposition - the wavelet
representation. Ieee T Pattern Anal 11:674–693
McDonald MG, Harbaugh AW (1988) A modular three-dimensional finitedifference ground-water flow model. USGPOMetropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E (1953) Equationof state calculations by fast computing machines. J Chem Phys 21:1087–1092
Meyer Y (1990) Wavelets, quadrature mirror filters and numerical imageprocessing. Lect Notes Math 1438:14–25
Montgomery DR, Manga M (2003) Streamflow and water well responses toearthquakes. Science 300:2047–2049
Pearson K (1901) On lines and planes of closest fit to systems of points in space.
Philos Mag 2:559–572
Petitta M, Mastrorillo L, Preziosi E, Banzato F, Barberio MD, Billi A, Cambi C, De
Luca G, Di Carlo G, Di Curzio D, Di Salvo C, Nanni T, Palpacelli S, Rusi S, Saroli
M, Tallini M, Tazioli A, Valigi D, Vivalda P, Doglioni C (2018) Water-table anddischarge changes associated with the 2016-2017 seismic sequence incentral Italy: hydrogeological data and a conceptual model for fractured
carbonate aquifers. Hydrogeol J 26:1009–1026
Roeloffs E (1996) Poroelastic techniques in the study of earthquake-relatedhydrologic phenomena. Adv Geophys 37(37):135–195
Roeloffs E, Quilty E (1997) Case 21: Water level and strain changes preceding andfollowing the August 4, 1985 Kettleman Hills, California, earthquake. Pure
Appl Geophys 149:21–60
Roeloffs EA (1988) Hydrologic precursors to earthquakes - a review. Pure Appl
Geophys 126:177–209
Rojstaczer S, Wolf S (1992) Permeability changes associated with large
earthquakes - an example from Loma-Prieta, California. Geology 20:211–214
Smith LC, Turcotte DL, Isacks BL (1998) Stream flow characterization and featuredetection using a discrete wavelet transform. Hydrol Process 12:233–249
Theis CV (1935) The relation between the lowering of the piezometric surface
and the rate and duration of discharge of a well using ground water storage.
Eos T Am Geophys Un 16:519–524
Ujiie K, Kimura G (2014) Earthquake faulting in subduction zones: insights fromfault rocks in accretionary prisms. Prog Earth Planetary Sci 1:7
Wakita H (1975) Water wells as possible indicators of tectonic strain. Science 189:
553–555
Wang CY, Cheng LH, Chin CV, Yu SB (2001) Coseismic hydrologic response of an
alluvial fan to the 1999 Chi-Chi earthquake, Taiwan. Geology 29:831–834
Wang C-Y, Liao X, Wang L-P, Wang C-H, Manga M (2016) Large earthquakes
create vertical permeability by breaching aquitards. Water Resour Res 52:
5923–5937
Wang C-Y, Wang C-H, Kuo C-H (2004a) Temporal change in groundwater level
following the 1999 (Mw = 7.5) Chi-Chi earthquake, Taiwan. Geofluids 4:210–
220
Wang C-Y, Wang C-H, Manga M (2004b) Coseismic release of water frommountains: Evidence from the 1999 (Mw = 7.5) Chi-Chi, Taiwan, earthquake.
Geology 32:769–772

Huang et al. Progress in Earth and Planetary Science

Authors’ contributions
Dr. Chien-Lin Huang proposed the topic and conceived and designed thestudy. Ph.D. student Chun-Hao Yao carried out the experimental study. Dr.
Chien-Lin Huang and Professor Nien-Sheng Hsu analyzed the data andhelped in their interpretation. Professor Nien-Sheng Hsu and Associate
Professor Gene J.-Y. You collaborated with the corresponding author in theconstruction of manuscript. All authors read and approved the final
manuscript.
Funding
This research was partially supported by the Ministry of Science andTechnology, Taiwan (Grant Nos. MOST 107-2625-M-002-005 and 108-2625-M002-007).
Availability of data and materials
The analyzed data was measured by the Water Resources Agency, Taiwan,
and acquired from the website https://gweb.wra.gov.tw/Hydroinfo/. Due tothe copyright and intellectual property rights of the data used in this study,
the readers can acquire the data through the above website or bycontacting the author for data requests.
Competing interests
The authors declare that they have no competing interests.
Received: 3 January 2020 Accepted: 25 June 2020

References
Adamowski J, Chan HF (2011) A wavelet neural network conjunction model forgroundwater level forecasting. J Hydrol 407:28–40
Bredehoeft JD (1967) Response of well-aquifer systems to earth tides. J Geophys
Res 72:3075
Cannas B, Fanni A, See L, Sias G (2006) Data preprocessing for river flow
forecasting using neural networks: wavelet transforms and data partitioning.
Phys Chem Earth 31:1164–1171
Chia Y, Chiu JJ, Chiang Y-H, Lee T-P, Wu Y-M, Horng M-J (2008) Implications ofcoseismic groundwater level changes observed at multiple-well monitoring
stations. Geophys J Int 172:293–301
Daubechies I (1988) Orthonormal bases of compactly supported wavelets.
Commun Pure Appl Math 41:909–996
Donoho DL (1995) De-noising by soft-thresholding. Ieee T Inform Theory 41:613–
627
Elkhoury JE, Brodsky EE, Agnew DC (2006) Seismic waves increase permeability.
Nature 441:1135–1138
Freeze RA, Cherry JA (1979) Groundwater. Prentice-Hall, Englewood Cliffs
Hartmann J, Levy JK (2006) The influence of seismotectonics on precursory
changes in groundwater composition for the 1995 Kobe earthquake, Japan.
Hydrogeol J 14:1307–1318
Hasegawa R, Yamaguchi A, Fukuchi R, Hamada Y, Ogawa N, Kitamura Y, Kimura
G, Ashi J, Ishikawa T (2019) Postseismic fluid discharge chemically recordedin altered pseudotachylyte discovered from an ancient megasplay fault: an
example from the Nobeoka Thrust in the Shimanto accretionary complex,
SW Japan. Prog Earth Planetary Sci 6
Jacob CE (1940) On the flow of water in an elastic artesian aquifer. Eos T Am
Geophys Un 21:574–586
Jang CS, Liu CW, Chia Y, Cheng LH, Chen YC (2008) Changes in hydrogeologicalproperties of the River Choushui alluvial fan aquifer due to the 1999 Chi-Chi
earthquake, Taiwan. Hydrogeol J 16:389–397
Kagabu M, Ide K, Hosono T, Nakagawa K, Shimada J (2020) Describing coseismic
groundwater level rise using tank model in volcanic aquifers, Kumamoto,southern Japan. J Hydrol 582
Katayama I, Kubo T, Sakuma H, Kawai K (2015) Can clay minerals account for the
behavior of non-asperity on the subducting plate interface? Prog Earth
Planetary Sci 2:30

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.