Identification of hydrogeological evolution using hydrogeology-seismology analysis of groundwater head fluctuation related to the 1999 MW = 7.5 Chi-Chi earthquake

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 the hydrogeological 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 average leaking/injecting rate, (2) distance between the acting position and monitoring well, (3) storage coefficient under effective 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 show that the overlying principal components PCs and low-level wavelet de-noising can filter additional sources/sinks, in which 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 pre-seismic 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 was increased to dissipate excess pore-water pressure and effective stress.


Introduction
The crustal stress (strain) caused by earthquakes is a natural driving force, causing the groundwater head to change in a short period in a large area subject to a consistent stressor (Wakita 1975). Hence, it provides a medium for understanding the hydrogeological evolution of the formation through earthquake-induced changes in groundwater head. The changes in groundwater head caused 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 elastic change caused by seismic wave transmission, in which the groundwater head returns to the original head after the seismic wave dissipates. The step-like changes in groundwater head are mainly affected by the inelastic change 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 the fact that the stratum is squeezed by the seismic stress, causing the water pressure of the large area to rise, so that the excess pore-water pressure may be dissipated to the surface. The co-seismic step-like change in groundwater head is usually accompanied by a post-seismic recovery phenomenon, which represents the pore-water pressure in the stratum stabilizing. This behavior is mainly controlled by the hydrogeological characteristics of the porous medium that can be identified from the dissipation 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), andHartmann andLevy (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) found that the regional distribution of co-seismic groundwater head rise and fall closely coincided with areas of contraction and dilatation expected from faulting. In research focusing on post-seismic hydrogeological changes, Rojstaczer and Wolf (1992), Montgomery and Manga (2003), and Petitta et al. (2018) found that earthquakes result in persistent or permanent changes in aquifer characteristics, hydrogeological parameters, and groundwater hydrology.
The 1999 Chi-Chi (M w 7.5) earthquake was the largest to hit Taiwan in the past century. A dense network of hydrological stations on the Chou-Shui River alluvial fan near the epicenter of the earthquake captured the groundwater head changes before, during, and after the earthquake. The proximity to a large earthquake and the dense network of hydrological stations provided a rare opportunity to examine the effect of an earthquake on groundwater. Several studies have discussed the abnormal groundwater head changes in response to this earthquake 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 and Kimura 2014), but few previous studies considered the earthquake-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 with excess pore-water pressure release and effective stress relief should be in evolution, so many parameters may increase the uncertainties and error of simulating the drawdown in groundwater head. Furthermore, previous studies have seldom used a stochastic approach to simulate all possible natural hydrogeological scenarios, in which the measured groundwater head may involve additional factors such as sources/sinks and boundary conditions, where the corresponding magnitude may be filtered by proper spatiotemporal-frequency analysis.
Accordingly, the purpose of this study is 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 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 strain caused 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 the identification of drawdown related to excess pore-water pressure 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.

Methodology
Procedures 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 groundwater head filtration; (3) establishment of the stochastic experimental optimization model; and (4) solution of the stochastic optimization model for identifying the evolving process. The flowchart of the methodology is shown in 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, this study firstly investigates three natural phenomena in hydrogeology-seismology: (1) vertical concussion hot zones 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 followed by relaxation from excess pore-water pressure. Three different signal decomposition approaches are used to evaluate the strain relaxation process, high-low-rank singular value decomposition (SVD), high-low-frequency wavelet transform decomposition, and high-low-level soft-thresholding de-noising. These approaches are used in conjunction with the measured earth-oceanic tidal head, barometric pressure, and other geophysical factors.
Step 3-1. Embed a post-seismic groundwater head drawdown 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 model by 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 is computed from a series of stochastic experimental calculated g α, b, i .
Step 4. Apply the simulated annealing algorithm to solve 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 and post-seismic evolution.
Global spatiotemporal investigation of hot zones pattern of change in co-post-seismic groundwater head and detection 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, in which the partial derivatives of groundwater head vs. time should be larger than 5 cm, i.e., j ∂h f ðt;κÞ ∂t |≥ 5 cm.
(2) Investigate the concussion change in co-postseismic groundwater head across multiple zones κ and aquifers f ( ∂h f ðt;κÞ ∂t ) by using quartile analysis and the spatial patterns of change in co-seismic groundwater head across multiple aquifers (i.e., ∂h f ðt;κÞ ∂t 0 ).
(3) Investigate the spatiotemporal pattern of postseismic uplift/downward partial derivatives of groundwater head vs. time during excess porewater pressure release periods t 1 and t 2 across multiple zones and aquifers (Εt½ ∂h f ðt;κÞ (4) Detect the hot zone patterns of obvious vertical permeability increase during excess pore-water pressure release evolving periods t 1 and t 2 across multiple zones and aquitards T (Εt½ ∂h f ðt;κÞ ), in which the aquitard T is coated by aquifer f and for the zone κ means that the groundwater head in one of the aquifers among f and f + 1 rises and that of 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 the step-like sudden rise/drop in co-seismic groundwater head followed by a release/recovery behavior is mainly caused by crustal strain followed by relaxation. However, the post-seismic-measured groundwater heads might be affected by factors (source/sink) such as rainfall R rain , artificial pumping P human , earth tidal T earth , oceanic tidal T ocean , barometric effect B baro , and other factors O (which include well loss, wellbore storage effects, and noise from automatic recording). These existing sources/ sinks, which are formulated as W in Eq. (1), will affect groundwater flux. The governing equation of the threedimensional groundwater flow can be expressed as (McDonald and Harbaugh 1988): where K xx , K yy , and K zz are the hydraulic conductivities along the x, y, and z directions, respectively 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.
Principal component analysis (PCA) Longuevergne et al. (2007) firstly applied PCA to decompose the groundwater head hydrograph into a series of intrinsic mode functions from low-rank to high-rank components 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 new variables (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 factor scoreC and the factor loading F as follows (Longuevergne et al. 2007): where the row vector of the n × κ f -order matrix U is an orthonormal left singular vector {u 1, · , f , …, u n, · , f }; The row vector of the κ f × κ f -order matrix V is an orthonormal right singular vector fv 1;Á; f ; …; v κ f ;Á; f g . The singular value σ j, f determines the eigenvalue λ j; f ¼ 1 n-1 σ 2 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 κ th f monitoring well's groundwater head and the j th principal component coefficient.
The mathematical meaning in Eq.
(2) can be interpreted as there being L independent factors causing the drawdown in groundwater head. The factor score (transpose) matrix C T ¼ ½c 1;Á; f ⋯ c n;Á; f stores n time serial observations for L = κ f factors in aquifer f, where c t;Á; f ¼ ðc t;1; f ; …; c t; j; f ; …; c t;L; f Þ T . Hence, Eq. (2) can be rewritten as: 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-wise factor score c t; 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 h post-sei-relax t;κ f ; f and the corre- can be generated as follows: 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 the approximation coefficient and reduces partially leveled detail coefficient to identify the groundwater head drawdown factors which are associated with the post-seismic pressure 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: 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 levels of the wavelet transform, as shown in Fig. 2b. For N ∈ ℕ, a Daubechies wavelet ψ(x) and scaling function φ(x) of class D-2N is defined as follows: where N is the number of Daubechies coefficients η k , and η 0 , …, η 2N − 1 ∈ ℝ are the constant filter coefficients satisfying the conditions: 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 associated with the post-seismic relaxation process. Suppose that there are n noisy samples of a function f de (t) and we observe data The most general de-noising problem for the noisy groundwater head signal is expressed as follows Donoho (1995): where f de (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 to suppress the signal noise X(t) and recover f de (t). Our goal is to estimatef de ¼ ½f de ð1Þ; …;f de ðtÞ; …;f de ðnÞ with small mean square error, as expressed below: whered ðν − lÞ k is the de-noised detail coefficient in level ν − l. The noise is concentrated in the high-frequency IMFs, so those are removed by using the fast discrete Daubechies wavelet transform with a soft threshold. The soft-thresholding rule is usually referred to as wavelet shrinkage since it reduces the high amplitude detail coefficients d ðν − lÞ k so that their absolute value is below a certain threshold level, denoted by λ(l, k), which can generally be a function of the decomposition/resolution level l and the index k. This study calculates the denoised wavelet coefficientsd ðν − lÞ k using soft-thresholding rule δ soft λ , which the associated threshold λ(l, k) is calculated according to Donoho (1995), as expressed below: where σ(l, k) is the standard deviation of the noise signal in level l. The last step is to compute wavelet recon- Formulations of co-post-seismic groundwater head recovery equation

Force point source analytical solution
According to Roeloffs (1996), considering a saturated porous medium in an infinite domain, assuming that the medium is homogeneous isotropic soil and is subjected to three concentrated forces produced by the earthquakeinduced crustal strain along the coordinate axis (F r , F θ , and F z ) 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 u r , u θ , and u z and the excess pore-water pressure P are the variables and the coordinate origin is placed on the force point source, the governing equation of the pore elasticity theory can be expressed in cylindrical coordinates: where r is the horizontal distance from the origin, z is the depth of the force point source, ε is the volume ν is the Poisson's ratio, γ w is the specific weight of water, n p 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 (F r = F Z = F), letting D t ' = R 2 /4Dt and C 0 = 1/8πGηS s , where S s = K/D and D = T/S is the hydraulic diffusivity which can be estimated from the transmissivity (T) and storage coefficient (S). After simplification, the point source force (F) can determine the change in groundwater head (Δh t ) at different positions: where R is the distance between the force point source and monitoring point ( R ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 þ z 2 p ), 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 in groundwater head was gradually removed as time increases.
However, during co-post-seismic pressure release, the hydrogeological properties in the porous media are evolving, so accurate estimations of G, η, S s , F, T, S, and coefficients of the error function in Eq. (19) for different periods are difficult and should be time-variant, and should not equal the pre-seismic values.
Reduction of parameter and uncertainty by using transformed pore elasticity-based equation In reality, when an earthquake or swarm earthquake occurs, the forces along the horizontal and vertical directions are usually nonuniform (F r ≠ F Z ≠ F) and the tangential force F θ cannot be ignored, because the source, depth, and scale of each earthquake are different. Besides, the seismic point source forces often cause formation fracture, which produces the post-seismic 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 groundwater head 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: where H 0 is the initial groundwater head, Qð F; r; zÞ is the 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 for excess 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 and formation compression, and W ðuÞ is the harmonic average well function. From Roeloffs (1996) to this section, although the transformed Theis equation cannot estimate the groundwater head 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 reliable occurrence 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 associated with the relaxation process using ensemble analysis scheme ID number α. This study designs the uncertainty/reliability (occurrence probability) of each generated 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 is more 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 g thresh surrounded minimum that is sent to evaluate the post-seismic release scenario through hydrogeology. This study uses the root mean square error (RMSE) as an objective function g α, b, i to calculate the fitting accuracy between the measured/filtered drawdown s 0 mea α ðtÞ and the calculated drawdown s 0 cal b;i ðt; Q; r; S; TÞ, as expressed below: where s 0 mea α ðtÞ is the drawdown associated with the pressure relaxation filtered using analysis scheme ID number α; s 0 cal b;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, respectively; and ρ thresh Á;Á;Á is the exceedance probability threshold corresponding to the objective function value threshold g thresh , in which the transformation between ρ thresh Á;Á;Á and g thresh is made by the ODC. From Eqs. (15) to (18), one can find that F is positively proportional to Q r ( F∝ Q r ). 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 the crustal strain-induced fault (leaking/injecting source) and monitoring well, as expressed below: and Q UB ρ ð F; r; zÞ are the lower and upper bounds of the harmonic 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 and monitoring well of the referenced and analyzed monitoring zone, respectively; and r LB ρ ð F; r; zÞ and r UB κ ð 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 of the 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 groundwater head 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 and probability 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 of two adjacent EPs such as S E κ ðP κ;ρ $ P κ;ρþ1 ; G; η; S s ; F; r; z ; SÞ and T E κ ð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 average, namely ð S κ ðP κ;ρ ;ÁÞþS κ ðP κ;ρþ1 ;ÁÞ 2 Þ and ð T κ ðP κ;ρ ;ÁÞþT κ ðP κ;ρþ1 ;ÁÞ 2 Þ. Then, to evaluate the mixed influence of pore-water pressure release, effective stress relief, formation compression, and leaking/injecting effect, the harmonic evolving parameters are calculated from the summation of the expected parameters of all different EPs, as derived below: where S HAE κ ðÁÞ and T HAE κ ðÁÞ are the post-seismic harmonic average storage coefficients under effective stress compression and the transmissivity for pressure release of 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.

Solving stochastic optimization model using simulated annealing algorithm
This study used a simulated annealing algorithm proposed by Kirkpatrick et al. (1983) to optimize/solve the minimum value of the objective function for different initial solutions i, boundary conditions b, and analytical schemes α. 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 a monitoring well, storage coefficient under effective stress relief-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 is shown in Fig. 3, and the optimization steps are described below: Step 1. Set initial solution i, boundary condition b, and analysis scheme α, and then set the initial temperature T 0 c ¼ 100 and initial decision variables ξ 0 ¼ ½Q; r; S; T 0 .
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 g 0 .
Step 2-1. Randomly generate a neighborhood decision variable ðξ 0 Þ τ ¼ ½Q 0 ; r 0 ; S 0 ; T 0 τ and input it into the postseismic groundwater head drawdown simulation model to calculate the analytical solution, then calculate the corresponding objective function value (g') τ for a searching iteration ID number τ.
Step 3. Examine whether ξ τ + 1 meets the stopping condition (i.e., the gradient of objective function g τ vs. iteration times ∂g ∂τ is lower than 1 × 10 −10 after 2000 iterations). If the stopping condition is satisfied, then output the optimized g value and the corresponding 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 τ is replaced by T τ 0 c for jumping out of local minima; otherwise, if it has not accepted, then set τ = τ + 1 and return to step 2-1.

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 from the Baguah Mountain and Dulliu Hill in the east to the Taiwan Straits in the west and from Wu River in the north to Peikang River in the south. The radius of the alluvial fan is approximately 40 km and the area is approximately 2079 km 2 . The Chou-Shui River alluvial fan can be divided into proximal, mid, and distal fan regions according to topography, geology, and strata, as shown in Fig. 4a.
The groundwater strata in the Chou-Shui River alluvial fan is approximately 330 m thick, which can be divided 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   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 Yanlin-2(YL2) monitoring wells in 2004, where only SC2 in aquifer 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), and Ganghou-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 were combined with the SC2 and GH3 for a total of nine wells, which the datasets for comparison and subsequent signal analyses of the pre-co-post-seismic change in groundwater head is illustrated in Fig. 5b.

Results and discussions
Investigated hot zones and pattern of change in co-postseismic groundwater head across multiple sub-fans and aquifers Initially, this study investigated the hot zones of intense change in co-post-seismic hourly groundwater head across multiple sub-fans through quartile analysis, in which the partial derivatives of groundwater head vs. time (Δt = 1 h) followed j ∂h f ðt;κÞ ∂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 and downstream middle fan are larger than those of the northern middle fan, respectively. This is because the northern upstream boundary comprises the Changhua Fault, and the Chiuchiungkeng Fault is only minorly connected with the 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 groundwater head 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 the stratum 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 and proximal fan shows both increase and decrease in coseismic groundwater head. This study investigated the quartile change patterns in co-post-seismic hourly groundwater head ( ∂h f ðt;κÞ ∂t ) across multiple sub-fans. Figure 7 indicates that the extreme values occur substantially across the proximal, middle, and distal fans, in which the rapid rise and 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 across multiple aquifers.
5 ) during 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 than that 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. Additionally, from Fig. 10, this study demonstrates that the partial derivative groundwater head in the continuously pressurized zones in Fig. 9 during post-seismic t = 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 and northern boundary in aquifer 4 continuously suffer crustal strain pressure.

Detected spatiotemporal hot zones of obvious vertical permeability increase
Hot zones of obvious vertical permeability increase across multiple aquitards and zones κ ( 3 5 ) are investigated by using the Kriging gridding method. Figure 11 shows that the condition Ε t ½ ∂h f ðt;κÞ ∂t Á ∂h f þ1 ðt;κÞ ∂t 5 < 0 is met in part of the middle and proximal fans in aquitards 1 and 2 and at the distal fan in aquitard 3 during post-seismic t = 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 it means 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-like sudden rise/drop in co-seismic groundwater head followed 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, R rain , 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, P human . 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-water pressure release/recovery behavior, and the other highrank decomposed drawdown may be related to other geophysical sources/sinks (earth tidal, barometric pressure, recharge) and noise. Therefore, this study used four combinations of PC1, PC1+PC2, PC1+PC3, and PC1+ 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 s 0 mea α in Eq. (25) for optimization.
Wavelet transformed pattern of co-post-seismic groundwater heads fluctuation One to three levels of wavelet transform were executed on the post-seismic-measured groundwater heads of SC2 and GH3 monitoring wells to calculate the approximation coefficients cA 1 , cA 2 , and cA 3 (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 s 0 mea α in Eq. (25) for optimization. The abovementioned figures demonstrate that the decomposed low-level drawdown hydrograph shows more refined drawdown than the high-level drawdown hydrograph, 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 performed on 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 the applied 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 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 a post-seismic pumping test. From Fig. 4, we see the connection between the epicenter and SC2 monitoring well and between the epicenter and GH3 well all pass through the Chelungpu Fault outside of the Chou-Shui River alluvial fan and Changhua Fault. This study uses the Changhua Fault as a reference to calculate 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 the well and fault, as expressed below: Moreover, according to the in situ hydrogeological soil composition of SC2 and GH3 in a confined aquifer, the corresponding boundary condition is set as expressed below: Furthermore, this study designs a series of experiments for 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 and effective stress delivered through hydrogeology in Si-Chou (2) After the stochastic experimental optimization, the evaluated evolving parameters (leaking/injecting rate Qð F; r; zÞ, harmonic average distance between the 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 Þ) for different initial conditions i, boundary conditions b, and analysis schemes α were collected. The stochastic optimized solutions for SðG; η; S s ; F; r; z; SÞ and T ð F; r; z; T Þ associated with the corresponding objective function value g α, b, i for SC2 is shown in Fig.  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 (Jang et al. 2008), and the evolving transmissivity increased from a pre-seismic value of 92.4 m 2 /h, post-seismic 27th hour value 98.6 m 2 /h to post-seismic value of 147.6 m 2 /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 shows that 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 associated with the excess pore-water pressure relaxation process, as shown in Fig. 16a, the optimized objective function value (RMSE 0.06327 m) was the lowest among the three filtration analysis methods. Ensemble optimized pattern of pore-water pressure and effective stress delivered through hydrogeology in Gang-Hou (3) The stochastic experimental optimized solutions for SðG ; η; S s ; F; r; z; SÞ and Tð F; r; z; T Þ associated with the corresponding 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 to release 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 27th hour value 0.000112, which is a decrease of 24.8%; and the evolving transmissivity increased from a pre-seismic value of 28.8 m 2 /h to post-seismic 27th hour value 120.7 m 2 /h, meaning that the GH3 zone also experienced large crustal compression and permeability was increased to dissipate excess pore-water pressure and effective stress in 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 from other filtration analysis methods. Furthermore, Fig. 14b-2 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, so that the S can be expressed as Eq. (30). Accordingly, it can be seen that when the earthquake compresses the aquifer, the soil porosity η and the coefficient of compressibility of the soil α decrease, so that S also decreases.
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 usually reflected by volume changes caused by the sliding (or dislocation) of aquifers composed of poroelastic media at seismic faults. If an aquifer is located at a crustal compression zone, its groundwater heads will rise; conversely, if the aquifer is located at a crustal expansion zone, 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 both showed a rise in co-seismic groundwater head and a fall in the post-seismic identified evolving storage coefficient Sð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 to damage 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 the occurrence 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 of zones 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 be understood 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 the tide data from the Mailiao (ML) station measured by the Central 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 cD 2 in Fig. 17b show the earth tidal fluctuation effect, T earth . Moreover, Jacob (1940) found that the changes in groundwater head and barometric pressure show opposite trends. The barometric data from Taichung (TC) station operated by the Central Meteorological Bureau is shown in Fig. 17d. Accordingly, the drawdown extracted from PC4+PC9 in Fig. 17a and cD 3 in Fig. 17b could be approximately regarded as the barometric effect, B. With regard to the drawdown extracted from PC8 in Fig. 17a and cD 1 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.

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-water pressure 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 spatiotemporal change in co-post-seismic groundwater head and vertical permeability 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 and multi-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 for excess 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 the stochastic optimization model. Moreover, this study establishes an indicator for evaluating the post-seismic hydrogeological harmonic average evolving parameters according to the Riemann integral and an objective duration curve.
At 1:47 a.m. on September 21, 1999, the M w 7.5 Chi-Chi earthquake struck Taiwan. This study applied the integrated 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 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 southern boundary 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 part of 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 geophysical factors 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 cD 2 could be approximately regarded as the earth tidal fluctuation effect. Moreover, the drawdown from PC4+PC9 and cD 3 could be understood as the barometric effect, and the drawdown from PC8 and cD 1 as white noise and other interference factors.
After the stochastic experimental optimization process, according to the Riemann integral and an ODC calculated occurrence probability for SC2, the storage coefficient was reduced from the pre-seismic on-site pumping test value 0.00107, post-seismic 27th hour 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 m 2 /h, post-seismic 27th hour value 98.6 m 2 /h, to post-seismic value 147.6 m 2 /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 m 2 / h to post-seismic 27th hour value of 120.7 m 2 /h. These results demonstrate that zones SC2 and GH3 suffered large crustal compression resulting in permeability increase to dissipate excess pore-water pressure and effective stress in the deep confined aquifer.