# Estimating the hydraulic parameters of a confined aquifer based on the response of groundwater levels to seismic Rayleigh waves

Estimating the hydraulic parameters of a confined aquifer based on the response of groundwater... Summary Groundwater flow models implemented to manage regional water resources require aquifer hydraulic parameters. Traditional methods for obtaining these parameters include laboratory experiments, field tests and model inversions, and each are potentially hindered by their unique limitations. Here, we propose a methodology for estimating hydraulic conductivity and storage coefficients using the spectral characteristics of the coseismic groundwater-level oscillations and seismic Rayleigh waves. The results from Well X10 are consistent with the variations and spectral characteristics of the water-level oscillations and seismic waves and present an estimated hydraulic conductivity of approximately 1 × 10−3 m s−1 and storativity of 15 × 10−6. The proposed methodology for estimating hydraulic parameters in confined aquifers is a practical and novel approach for groundwater management and seismic precursor anomaly analyses. Permeability and Porosity Hydrogeophysics, Hydrology, Fourier analysis, Time series analysis 1 INTRODUCTION Understanding hydrological properties is crucial for many applications, such as groundwater resource management, subsurface solute transport assessments, hydraulic fracturing and mining. Transmissivity (T) and storativity (S) are the two key parameters that must be accurately determined in many of these applications. The methods for estimating T and S include laboratory experiments, field tests and water-level responses in wells caused by natural perturbations. Laboratory methods analyse the physical properties (such as the porosity, grain size and compression coefficient) of the rock or soil samples (Chapuis 1992; Ma et al.2014). Field tests are conducted by pumping a well and then observing water-level responses in the pumped well or in nearby observation wells. Data from these tests are interpreted using a range of analytic solutions, including the Theis, corrected Theis and Hantush models (Theis 1935; Jacob 1940; Cooper & Jacob 1946). The complex geometry, boundary conditions and heterogeneity parameterizations have been incorporated into numerous analytical pump test solutions (Dawson & Istok 1991; Jiao & Rushton 1995; Sethi 2011). Water-level responses to natural perturbations include barometric-pressure fluctuations, surface loading and earth tides (Hsieh et al.1987; Allègre et al.2016). Of these three methods, laboratory methods may provide inaccurate parameter values because of the natural heterogeneity and artefacts associated with the specific testing protocol. Although field testing is common, it can be expensive, and continuous observations of monitoring wells may not be feasible. The response to natural perturbations is becoming more widely used because real-time monitoring systems are increasingly adopted and appropriate methods for determining hydraulic parameters have been developed (Xue et al.2013; Wang et al. 2016; Shi & Wang 2016), especially the responses of well-water levels to earth tides to determine the aquifer hydraulic conductivity. However, this method is hindered because it requires a clear tidal signal in the observation wells. Water-level responses to earthquakes represent another method for estimating aquifer properties (Shih 2009; Sun et al.2015; Wang et al.2016). The main coseismic groundwater-level response types are step change and oscillation change (Shi et al.2015). Coseismic step changes have been used to estimate permeability (Wang et al.2016), and oscillation changes have been used to estimate specific storage (Shih 2009; Sun et al.2015). However, to the best of our knowledge, a method of estimating aquifer hydraulic properties from coseismic responses has not been proposed. Therefore, we present a method that compares the correlation between the spectral characteristics of coseismic water-level oscillation and seismic Rayleigh waves, which we then use to identify the hydraulic properties of a well in a confined aquifer. 2 BACKGROUND AND OBSERVATIONS Well X10 is located to the south of Urumqi City in the Xinjiang Autonomous Region in western China at 43.70° N and 87.62° E and an elevation of 1056 m amsl. The well hole was drilled in September 1980. Well X10 is at the intersection of the Liushugou-Hongyanchi thrust fault and its derived fault (Fig. 1). Groundwater discharges through springs at this location provide stable year-round yields of 20–30 l s−1. Figure 1. View largeDownload slide Geological setting, location of Well X10 (blue triangle) and epicentres of the 10 earthquakes (beach balls) that occurred between November and December 2016 Figure 1. View largeDownload slide Geological setting, location of Well X10 (blue triangle) and epicentres of the 10 earthquakes (beach balls) that occurred between November and December 2016 As shown in Fig. 2, the total depth of Well X10 is 28 m, the diameter is 146 mm at depths to 9.2 m, the diameter is 130 mm at depths from 9.2 to 13.04 m, and the diameter is 110 mm at depths from 13.04 m to the bottom. Core samples collected between 16 and 24 m showed extensive fracturing; therefore, a screen filter was installed in this zone. The aquifer in the Permian Hongyanchi group is 8 m thick, consists of sandstone and gravel and has a static water level of approximately 1 m below ground surface. Figure 2. View largeDownload slide Aquifer lithology and wellbore structure of Well X10. Figure 2. View largeDownload slide Aquifer lithology and wellbore structure of Well X10. Starting in October 2016, the water level was monitored using an SWY-II type digital water-level metre, which was developed and produced by the Institute of Crustal Dynamics, China Earthquake Administration. The water-level sensor has a range of 0–50 m, resolution of 1 mm, sampling interval of 1 s, precision of 0.02 per cent and accuracy of 0.05 per cent F.S. Ten earthquakes during November and December 2016 caused the water level to fluctuate in Well X10 (Fig. 1). The water-level changes in a well induced by far-field earthquakes were strongly correlated to the dynamic strains. Wang (2007) introduced a physical quantity, the seismic energy density e (J m−3), which is a function of the earthquake magnitude (M) and epicentral distance, r (km): log r = 0.48M − 0.33 log e – 1.4. A comparison with other well sites and other earthquake-related hydrological phenomena, we plotted the earthquake events that induced water-level changes using the magnitude versus epicentral distance to Well X10 (Fig. 3). A global data set of water-level responses compiled by Wang & Manga (2010) were also plotted for comparison. The seismic energy densities were significantly lower in Well X10 than in the global data set, which were mainly within 10−4–10−7 J m−3, indicating that Well X10 is highly sensitive to seismic waves. Figure 3. View largeDownload slide Seismic energy density as a function of epicentral distance and earthquake magnitude. Lines of constant seismic energy density are given in J m−3, and the red squares and grey circles correspond to Well X10 and a global data set compiled by Wang & Manga (2010), respectively. Figure 3. View largeDownload slide Seismic energy density as a function of epicentral distance and earthquake magnitude. Lines of constant seismic energy density are given in J m−3, and the red squares and grey circles correspond to Well X10 and a global data set compiled by Wang & Manga (2010), respectively. The fluctuations in Well X10 are compared with vertical displacement waveforms observed by the SLJ-100 seismic recoding sensor at station SMG, which is located 10 km northeast of Well X10. Fig. 4 shows a comparison of the coseismic responses of the water levels in Well X10 with the vertical displacement waveforms at station SMG caused by the M 7.8 Solomon Islands earthquake (10.665° S, 161.335° E), which occurred at 01:38 China Standard Time (CST) on 2016 December 9. Figure 4. View largeDownload slide Water level changes in the Well X10 (a) and vertical displacement at SMG (b) due to the 2016 December 9 M7.8 Solomon Islands earthquake and time–frequency images of water level (c) and displacement (d). The colour bar shows the amplitude of water level or displacement. Figure 4. View largeDownload slide Water level changes in the Well X10 (a) and vertical displacement at SMG (b) due to the 2016 December 9 M7.8 Solomon Islands earthquake and time–frequency images of water level (c) and displacement (d). The colour bar shows the amplitude of water level or displacement. To better illustrate the similarities between the two waveforms, the S-transform method (Stockwell et al. 1996) was used to calculate the time–frequency characteristics as presented in Fig. 4(c) (water-level oscillation in Well X10) and Fig. 4(d) (seismic waves at the SMG station). Time is on the horizontal axis, frequency is on the vertical axis, and the colours represent response amplification. Note that the similarities in the shapes and frequencies of the water-level oscillations in Well X10 and the vertical component of the seismic waves at the SMG station were caused by the M7.8 Solomon Islands earthquake. The S-waves and surface waves were also clearly detected. A comparison of the time–frequency characteristics in Figs 4(c) and (d) show that the water and seismic waves also have a wide frequency range of 0.03–0.06 Hz (15–35 s) and last approximately 20 min. For the surface wave, the low frequency and long period signals were first observed and followed by high frequency and short period signals. The comparison of the amplifications shows that the water and seismic waves reached a range of values at a period of approximately 20–40 s, which was the main period of the surface wave (Rayleigh wave). 3 ESTIMATION OF HYDRAULIC PARAMETERS The analysis presented above highlights the similarity between the shapes and frequencies of the water-level oscillations and the vertical component seismic waves, especially for the surface wave period. Cooper et al. (1965) investigated the responses of aquifer systems to vertical surface movements and stated that the water-level responses were induced by both the vertical surface movements and pore-pressure fluctuations in the aquifer (Fig. 5), and their results also demonstrated that the response amplifications varied with the signals at different periods and were closely correlated with the hydraulic parameters (e.g. transmissivity and storativity). Figure 5. View largeDownload slide Idealized representation of an open well in which water level fluctuations are caused by the oscillation of artesian pressure and vertical motion of land surface. Where p = ρg(H + z) + p0 sin (ωt − η) represents the pressure fluctuation in aquifer with a frequency ω and phase lag η; h = h0 sin (ωt) means the water-level fluctuation caused by pressure fluctuation; and s = s0 sin (ωt − ζ) represents the fluctuation caused by surface land motion with a phase lag ζ. Figure 5. View largeDownload slide Idealized representation of an open well in which water level fluctuations are caused by the oscillation of artesian pressure and vertical motion of land surface. Where p = ρg(H + z) + p0 sin (ωt − η) represents the pressure fluctuation in aquifer with a frequency ω and phase lag η; h = h0 sin (ωt) means the water-level fluctuation caused by pressure fluctuation; and s = s0 sin (ωt − ζ) represents the fluctuation caused by surface land motion with a phase lag ζ. As a result, hydraulic parameters can be estimated by analysing the water levels during the surface wave period and the magnitudes of the ground motions (Sun et al.2015). Shih (2009) used the spectral analysis method to investigate the water response to seismic surface waves and estimated the storativity from spectral correlations of the water levels and the vertical displacement of the earth's surface. These methods are used to estimate aquifer hydraulic parameters by analysing the characteristics of the water-level oscillations in Well X10 in response to 10 earthquakes from around the world. 3.1 Storativity and specific storage The storativity (S, dimensionless) is defined as S = d·Ss, where d is the thickness of aquifer and Ss is defined as the volume of water released from a unit volume of porous aquifer when there is a unit decline in the head (Freeze & Cherry 1979), in the dimensions of length−1. The amplitude of the response of water levels to ground motions is closely correlated with aquifer thickness, storativity and surface wave number (Shih 2009):   $$h(t) = \frac{{0.918d{k_x}}}{S}w(t),$$ (1)where h(t) is the water level (dimensions of length), w(t) is the ground displacement (dimensions of length), and kx is the wave number on the surface. The amplitudes of the responses to signals at different periods were calculated using Fourier transforms (Duhamel & Vetterli 1990):   $$\begin{array}{@{}l@{}} {Z_h}(\omega ) = \int_{{ - \infty }}^{\infty }{{h(t){{\rm{e}}^{ - 2\pi i\omega }}^t{\rm{d}}t}}\\ {Z_w}(\omega ) = \int_{{ - \infty }}^{\infty }{{w(t){{\rm{e}}^{ - 2\pi i\omega }}^t{\rm{d}}t}}, \end{array}$$ (2)where Zh(ω) and Zw(ω) are the Fourier transforms of h(t) and w(t), respectively, ω is the frequency (dimensions of angular per time), t is time, and i = $$\sqrt { - 1} \$$. The time-series and frequency spectra of the water level and vertical displacement induced by the M7.8 Solomon Islands earthquake are plotted in Fig. 6. Compared with the waves at other periods, the response amplitude at the 20–30 s period was relatively large. Figure 6. View largeDownload slide Time-series and its corresponding spectral of water level and vertical displacement at the surface wave period. Figure 6. View largeDownload slide Time-series and its corresponding spectral of water level and vertical displacement at the surface wave period. The auto-spectral density (PSD) is used to detect the strong signal in time-series, whereas the cross-spectral density (CSD) and coherence measure the intensity of specific target signals between two time-series. The 95 per cent confidence interval of the PSD and non-zero coherence level for the CSD are also evaluated to identify significant components in the frequency domain (Shih et al.2013). Aquifer storativity was calculated from the spectral correlations of the water level and vertical displacement as follows (Shih 2009):   $$S = \frac{{1.836\pi d}}{\lambda }{\left( {\frac{{{S_{hh}}}}{{{S_{ww}}}}} \right)^{ - 1/2}}$$ (3)  $$S = \frac{{1.836\pi d}}{\lambda }{\left( {\frac{{{S_{hh}}}}{{{S_{hw}}}}} \right)^{ - 1}}$$ (4)  $$S = \frac{{1.836\pi d}}{\lambda }{\left( {\frac{{{S_{wh}}}}{{{S_{ww}}}}} \right)^{-1}}$$ (5)where λ is the wave length (dimensions of length), which can be determined from the seismic velocity and period; Shh and Sww are the PSD of the water level h(t) and the vertical displacement of ground motion w(t), respectively; and Swh and Shw are the CSDs of h(t) and w(t), respectively. PSD and CSD can be found using a finite length (L) of the observation data (Bendat & Piersol 2000):   \begin{eqnarray} {S_{hh}}(\omega) &=& \mathop {\lim }\limits_{L \to \infty } \frac{2}{L}E\left[ {|{Z_h}(\omega ,L){|^2}} \right]\nonumber\\ {S_{ww}}(\omega) &=& \mathop {\lim }\limits_{L \to \infty } \frac{2}{L}E\left[ {|{Z_w}(\omega ,L){|^2}} \right]\\ {S_{hw}}(\omega) &=& \mathop {\lim }\limits_{L \to \infty } \frac{2}{L}E\left[ {{Z_w}(\omega ,L) \cdot {Z_w}(\omega ,{L_r})} \right],\nonumber \end{eqnarray} (6)where the expected value E is the average value in a sampling period. The PSD and CSD of the water level and vertical displacement induced by the M7.8 Solomon Islands earthquake are plotted in Figs 7(a) and (b), respectively. The PSD with the 95 per cent confidence interval has lower and the upper extremes of 0.7364 and 3.8797 cm2 Hz−1, respectively (Fig. 7a). The non-zero coherence significance level is 0.68 (Fig. 7b). The plots indicate that the spectral densities of the water level and vertical displacement were largest and highly correlated in the 0.03–0.04 frequency range and could be used to estimate the storativity (Shih et al.2013). Figure 7. View largeDownload slide Autospectral density of water level and vertical displacement (a), and cross-spectral density (b), the short dash lines indicate the 95 per cent non-zero coherence level. Figure 7. View largeDownload slide Autospectral density of water level and vertical displacement (a), and cross-spectral density (b), the short dash lines indicate the 95 per cent non-zero coherence level. The water bearing layer is at depths of 16–24 m in Well X10, and its thickness is d = 8 m. The wave velocity v in the formula for the wavelength calculation λ = v/ω is obtained by dividing the distance to the earthquake epicentre by the difference between the arrival time of surface wave and the start time of water-level oscillation. The frequency ω was the highest point of the PSD in Fig. 7(a). According to eqs (3)—(6) and the spectral correlations of the water level and vertical displacement induced by the M7.8 Solomon Islands earthquake, the average estimated storativity in the Well X10 was 19.2 × 10−6. 3.2 Aquifer porosity Porosity (n) is the pore volume per unit aquifer volume. Jacob (1940) derived an equation that relates specific storage to porosity:   $$n = \frac{{{B_e}{S_s}}}{{\beta \rho }},$$ (7) where Be is the barometric efficiency (dimensionless), β = 4.6 × 10−10 m2 N−1 is the water compressibility and ρ = 103 kg m−3 is the fluid density. The barometric efficiency (Be) was estimated using Clark's method (Clark 1967) and the regression deconvolution method (Rasmussen & Crawford 1997). As shown in Fig. 8(a), the water level in Well X10 was influenced by the barometric pressure. We calculated the barometric efficiency based on the changes (differential values) in water level and barometric pressure and obtained the aquifer response from both the differential values and detrended values. The results show that the barometric efficiency and the aquifer response decrease gradually as the time interval increases when using the differential values, and the aquifer response is stable when using the detrended values (Fig. 8b). The water-level response to barometric pressure based on the differential values could reveal the short-term response, and the detrended values could reveal a long-term response. Therefore, we concluded that the barometric efficiency of Well X10 varies with the barometric-pressure period. Figure 8. View largeDownload slide Comparison of water level and barometric pressure (a) and the barometric efficiency/response corresponding to different time (b) in Well X10 between 2016 November 1 and December 31. Be is the barometric efficiency with different time interval, and Br is the barometric response with different time lag. The segment of Br corresponding to the time between 1 and 500 min was calculated by minute values, and the segment of Br corresponding to the time between 60 and 10 000 min was calculated by hourly values. Figure 8. View largeDownload slide Comparison of water level and barometric pressure (a) and the barometric efficiency/response corresponding to different time (b) in Well X10 between 2016 November 1 and December 31. Be is the barometric efficiency with different time interval, and Br is the barometric response with different time lag. The segment of Br corresponding to the time between 1 and 500 min was calculated by minute values, and the segment of Br corresponding to the time between 60 and 10 000 min was calculated by hourly values. The time intervals used in the method were selected such that the barometric-pressure-independent water-level change was as small as possible compared with the barometric-pressure change. Focusing the time-interval selection around the correlation coefficient maximizes the value of the barometric-pressure change while minimizing the barometric-pressure-independent water-level change. As shown in Fig. 8(b), the barometric efficiency value was 0.5036, which was the mean of the barometric efficiencies with a correlation coefficient larger than 0.75. In combination with the specific storage Ss, the porosity of the aquifer can be estimated. 3.3 Transmissivity and hydraulic conductivity The response amplitude of the water-level oscillations induced by seismic waves is closely related to the aquifer parameters, borehole structure and seismic wave periods. In addition to causing vertical fluctuations in the surface medium, seismic waves can also cause pore-pressure fluctuations in the aquifer. As shown in Fig. 5, the oscillations of the well-water level are affected by both the vertical motions of the earth's surface and the pore-pressure changes, which are in turn amplified by the water levels. Therefore, the response amplitude of the water-level oscillations is affected by both the pore-pressure fluctuations in well-aquifer systems and vertical motions of the earth's surface. In a confined well-aquifer system, the response of well-water levels to the pore pressure and vertical surface displacement can be represented by an amplification factor (or amplitude ratio) (Cooper et al.1965):   \begin{eqnarray} \left. \begin{array}{@{}l@{}} A = {\left[{\left(1 - \displaystyle\frac{{\pi r_w^2}}{{T\tau }}Kei{\alpha _w} - \frac{{4{\pi ^2}{H_e}}}{{{\tau ^2}g}} \right)^2} + {\left(\frac{{\pi r_w^2}}{{T\tau }}Ker{\alpha _w}\right)^2}\right]^{ - 1/2}}\\ A' = A \cdot \displaystyle\frac{{4{\pi ^2}{H_e}}}{{{\tau ^2}g}},{\alpha _w} = {r_w}{(\omega S/T)^{1/2}} \end{array} \right\}\nonumber\\ \end{eqnarray} (8)where A (dimensionless) is the water-level amplification factor of the aquifer's pore-pressure fluctuations, A’ (dimensionless) is the water-level amplification factor of the vertical surface displacements, rw is the wellbore radius(dimensions of length), He is the effective height of water column (dimensions of length), T is the transmissivity (dimensions of square length per time), S is the storativity, τ is the period of seismic wave (dimensions of time), ω is the angular frequency, g is the gravitational acceleration, and Kei et Ker are the imaginary and real parts of the zero-order Kelvin functions. If the height of the water-level oscillations induced by fluctuations of the aquifer's pore pressure is h and the height of the water-level oscillations induced by the vertical surface displacement is s, then the total height of the water-level oscillations is x = h + s. If R is the ratio of the height induced by pore pressure to the height induced by the vertical displacement and M is the amplitude ratio of the total height of the water-level to the vertical displacement, then   $$R = \frac{h}{s} = \frac{{2.7{E_w}}}{{\rho gnv\tau }},M = \frac{x}{s} = \left(\frac{x}{h}\right) \cdot \left(\frac{h}{s}\right) = A \cdot R.$$ (9) Water at 18° centigrade has a bulk modulus of elasticity of Ew = 2.2 × 109 N m−2 and a density ρ of 103 kg m−3 (Cooper et al.1965); the gravitational acceleration is g = 9.81 m s−2; v and τ are the wave velocity and wave period, respectively; and n is the aquifer porosity. R and M were calculated using eqs (8) and (9) and the observed amplitudes of the water level and vertical displacement. If the ratio of the amplification factors A/A’ = R’, the ratio of the fluctuations R has an inverse relationship with the seismic wave period τ, although R’ has a positive quadratic relationship with τ. The calculation of R requires the wave velocity v, which can be obtained from the distance from the earthquake epicentre divided by the time difference between the arrival time of surface wave and the start time of the water-level oscillations. The period τ can be determined from the time–frequency diagram in Fig. 6. The effective water column height He required by the R’ calculation is He = H + 3d/8 (Cooper et al.1965), where H is the well's water column height above the upper confining bed. Because the upper confining bed in Well X10 is at a depth of 16 m and the depth to the water table is 1 m, the thickness of the well's water column is 15 m; therefore, He = H + 3d/8≈18 m. The M values at every period τ can be obtained statistically from the time–frequency values, and the respective values of A and A’ can be calculated with their R values (Fig. 9). The values represented by green circles in Fig. 9 were calculated using the spectral amplitude ratio of water level and displacement shown in Fig. 6. The magenta solid lines in Fig. 9 were the theoretical curves that fit the spectral amplitude ratio. To identify the optimal hydraulic conductivity (K) that can best fit the observed data, one thousand groups of theoretical values (A and A’) determined by eq. (8) were calculated with different values of K (from 0.001 to 1 cm s−1), and the standard deviations of each group were obtained. Then, the optimal K was selected from the group with the smallest standard deviation. Figure 9. View largeDownload slide Observed (circles) and estimated (black solid line) amplification factors based on the water-level oscillations in the Well X10: (a) due to pore pressure fluctuations; (b) due to vertical displacement. Figure 9. View largeDownload slide Observed (circles) and estimated (black solid line) amplification factors based on the water-level oscillations in the Well X10: (a) due to pore pressure fluctuations; (b) due to vertical displacement. A wellbore radius of rw = 54 mm was used in the calculation. Other parameters, such as storativity (S), porosity (n), wave velocity (v), etc., were all obtained using eqs (3) to (7). From the M7.8 Solomon Islands earthquake data, we obtained S = 19.2 × 10−6, n = 0.25, and hydraulic conductivity K = 0.375 cm s−1. 4 RESULTS AND DISCUSSIONS We can also acquire aquifer parameters based on the coseismic responses induced by nine other earthquakes (Fig. 1) using the methods described above. The results are shown in Table 1. Table 1. Estimated hydraulic parameters of Well X10 aquifer. Quantity  Unit  2016/12/25  2016/12/20  2016/12/17  2016/12/9  2016/12/8  2016/12/7  2016/11/25  2016/11/25  2016/11/22  2016/11/13  Magnitude  Mw  7.6  5.5  7.9  7.8  6.0  6.6  6.6  6.9  6.9  7.8  Distance  km  18,537  775  8451  9543  113  4348  1248  13,818  4520  12,798  Amplitude of WL  mm  15.0  5.0  4.0  27.5  30.0  3.5  4.5  24.5  18.5  15.0  Amplitude of VD  mm  0.582  0.057  0.109  1.046  0.430  0.094  0.168  0.604  0.488  0.582  Period  s  23  8  21  28  10  17  26  32  28  18  Wave velocity  m s−1  3461  3598  3732  3579  4302  3641  3821  3917  3725  3072  Wave Length  m  80 542  28 780  79 613  101 796  44 051  62 134  97 828  125 333  105 955  56 166  Shh  10−5  1923 983  6699  84 874  5495 827  324 297  98 186  454 044  158 855  1067 488  1944 934  Sww  10−5  1394  1  44  9813  186  26  262  270  1484  784  Shw  10−5  48,108  58  1853  223,039  6729  1478  10,003  6181  39,271  36,008  Storativitya  10−6  15.42  22.39  13.25  19.15  25.09  12.15  11.32  15.18  16.24  16.49  Storativityb  10−6  14.33  13.96  12.66  18.40  21.74  11.18  10.39  14.33  16.02  15.21  Storativityc  10−6  16.60  35.89  13.88  19.94  28.96  13.21  12.34  16.08  16.46  17.88  Storativityd  10−6  15.45  24.08  13.26  19.16  25.26  12.18  11.35  15.19  16.24  16.53  Specific Storage  10−6 m−1  1.93  3.01  1.66  2.40  3.16  1.52  1.42  1.90  2.03  2.07  Porositye  —  0.20  0.32  0.18  0.25  0.33  0.16  0.15  0.20  0.22  0.22  Hydraulic  cm s−1  0.05  0.10  0.09  0.38  0.09  0.09  0.07  0.10  0.12  0.08  conductivityf                        Transmissivity  m2 d−1  325  691  588  2592  622  608  449  677  802  518  Hydraulic  m2 s−1  243  332  513  1565  285  578  458  516  571  363  diffusivity  Quantity  Unit  2016/12/25  2016/12/20  2016/12/17  2016/12/9  2016/12/8  2016/12/7  2016/11/25  2016/11/25  2016/11/22  2016/11/13  Magnitude  Mw  7.6  5.5  7.9  7.8  6.0  6.6  6.6  6.9  6.9  7.8  Distance  km  18,537  775  8451  9543  113  4348  1248  13,818  4520  12,798  Amplitude of WL  mm  15.0  5.0  4.0  27.5  30.0  3.5  4.5  24.5  18.5  15.0  Amplitude of VD  mm  0.582  0.057  0.109  1.046  0.430  0.094  0.168  0.604  0.488  0.582  Period  s  23  8  21  28  10  17  26  32  28  18  Wave velocity  m s−1  3461  3598  3732  3579  4302  3641  3821  3917  3725  3072  Wave Length  m  80 542  28 780  79 613  101 796  44 051  62 134  97 828  125 333  105 955  56 166  Shh  10−5  1923 983  6699  84 874  5495 827  324 297  98 186  454 044  158 855  1067 488  1944 934  Sww  10−5  1394  1  44  9813  186  26  262  270  1484  784  Shw  10−5  48,108  58  1853  223,039  6729  1478  10,003  6181  39,271  36,008  Storativitya  10−6  15.42  22.39  13.25  19.15  25.09  12.15  11.32  15.18  16.24  16.49  Storativityb  10−6  14.33  13.96  12.66  18.40  21.74  11.18  10.39  14.33  16.02  15.21  Storativityc  10−6  16.60  35.89  13.88  19.94  28.96  13.21  12.34  16.08  16.46  17.88  Storativityd  10−6  15.45  24.08  13.26  19.16  25.26  12.18  11.35  15.19  16.24  16.53  Specific Storage  10−6 m−1  1.93  3.01  1.66  2.40  3.16  1.52  1.42  1.90  2.03  2.07  Porositye  —  0.20  0.32  0.18  0.25  0.33  0.16  0.15  0.20  0.22  0.22  Hydraulic  cm s−1  0.05  0.10  0.09  0.38  0.09  0.09  0.07  0.10  0.12  0.08  conductivityf                        Transmissivity  m2 d−1  325  691  588  2592  622  608  449  677  802  518  Hydraulic  m2 s−1  243  332  513  1565  285  578  458  516  571  363  diffusivity  aFrom Shh/Sww. bFrom Shh/Shw. cFrom Shw/Sww. daverage. eFrom Clark's method and eq. (7). fFrom Cooper's model. Notes: WL: water level; VD: vertical displacement. View Large Although the magnitudes and epicentres of the 10 earthquakes differ (the epicentral distances vary from 8451 to 18 537 km, which is almost half of the earth's perimeter), the estimated hydraulic parameters are similar. The porosity ranges from 0.15 to 0.33, the storativity ranges from 11.4 × 10−6 to 25.3 × 10−6 (corresponding to a specific storage Ss = S/d of 1.42–3.16 × 10−6 m−1), and the hydraulic conductivity ranges from 0.05 to 0.12 cm s−1 (corresponding to a transmissivity T = K·d of 300–800 m2 d−1). The hydraulic conductivity following the M7.8 earthquake on 2016 December 9 is 0.375 cm s−1, which is much larger than the other nine earthquakes. 4.1 Aquifer parameter changes induced by earthquakes Many researchers have noted that the aquifer permeability can be enhanced by seismic waves (Brodsky et al.2003; Elkhoury et al.2006; Manga et al.2012; Wang et al.2016), and this has also been proven by rock physics experiments (Elkhoury et al.2011; Candela et al.2014). Aquifer permeability enhancement always occurs with the coseismic water level abrupt/gradual response (Brodsky et al.2003; Sun et al.2015). The water-level changes in Well X10 were induced by the earthquakes from November—December 2016; however, nine of the ten showed oscillations but no subsequent step changes, indicating that the aquifer permeability did not change following the seismic waves (Shi et al.2015). A step-rise change occurred in the corrected water level following the M6.0 earthquake on 2016 December 8 (Fig. 8a), although the hydraulic conductivity did not increase immediately and showed an increasing water level after the earthquake. As a result, the estimated hydraulic conductivity during the M6.0 earthquake on 2016 December 8 was similar to the other earthquakes, and it showed substantial enhancement during the M7.8 earthquake on 2016 December 9, which was significantly higher than the other nine earthquakes. We speculated that the highest hydraulic conductivity estimated from the M7.8 earthquake data may be caused by the enhancement of hydraulic conductivity in the Well X10 aquifer after the M6.0 earthquake on 2016 December 8. Such enhancement, however, was caused by seismic wave shaking and new fractures induced by the stress/stain adjustment after the M6.0 earthquake (113 km from Well X10). The distribution of the coseismic strain field shows that Well X10 was in the compressive area of M6.0 earthquake and subjected to a coseismic strain of approximately 6 × 10−9 (Fig. 10), and the earthquake also has a higher seismic energy density than the other nine earthquakes (Fig. 3). As shown in Fig. 1a, Well X10 was drilled into a fault zone, which is especially sensitive to stress–strain variations. Figure 10. View largeDownload slide Coseismic strain field of the 2016 December 9 M7.8 Solomon Islands earthquake. Figure 10. View largeDownload slide Coseismic strain field of the 2016 December 9 M7.8 Solomon Islands earthquake. The coseismic strain was too small to increase the hydraulic conductivity, but the higher seismic energy density or the stress-strain adjustment after the M6.0 earthquake caused by the closer distance has the potential to create new fractures in the fault zone, which subsequently increases the porosity and hydraulic conductivity. As shown in Fig. 8(a), the step-rise recovery in the corrected water level (corresponding to the parameter changes) lasts for approximately 10 d and did not recover when the M7.8 earthquake occurred (12 hr after the M6.0 earthquake). 4.2 Model sensitivity to the estimated parameters The parameter estimations show that the storativity and hydraulic conductivity were less controlled by the earthquake magnitude and epicentral distance (Fig. 11, seismic energy density is closely related to the earthquake magnitude and epicentral distance) while the storativity was affected by the main seismic wave period. The storativity estimated from two earthquakes with smaller periods, the M6.0 earthquake on 2016 December 8 and M5.5 earthquake on 2016 December 20, was significantly greater than that for the other earthquakes, which may be related to the closeness between the significant seismic wave period (the corresponding period of the highest CSD at 10 s and 8 s of M6.0 and M5.5, respectively) and the resonant period of Well X10 $$(2\pi \sqrt {{H_e}/g} = 8.5\,{\rm{s}})$$. The amplification factor of the well-aquifer system will increase when the seismic wave period is closer to the well resonant period, which will increase the ratio of the power spectrum, Shh/Sww, Shh/Sww, and average storativity. However, the seismic wave period would not affect the hydraulic conductivity estimate because data from the period ranging from 5 to 50 s were all collected to estimate the hydraulic conductivity (Fig. 9), which was different from the estimation of storativity (only the significant period was selected). In addition, the highest estimated hydraulic conductivity in Fig. 11b on 2016 December 9 was caused by the enhancement of aquifer permeability after the M6.0 earthquake on 2016 December 8. Figure 11. View largeDownload slide The relation between seismic energy density and storativity (a), hydraulic conductivity (b), in which the colour circle shows the period of seismic wave. Figure 11. View largeDownload slide The relation between seismic energy density and storativity (a), hydraulic conductivity (b), in which the colour circle shows the period of seismic wave. Eqs (3)–(6) show that the storativity S was mainly affected by the wave length λ, aquifer thickness d and finite length L, where λ and d were determined by the observed data and L needs a manual setting. As shown in Fig. 12(a), different L values result in different curves of CSD, and the main difference is the smoothness of those curves. Although the smoothness was affected by the manual setting of L, the estimated results of porosity n and storativity S did not change significantly. Figure 12. View largeDownload slide Model sensitivity to the estimated parameters: (a) curve smoothness of CSD caused by different length (L) of data, and the corresponding porosity (n) and storativity (S); (b) estimated values of hydraulic conductivity (K) by inputting different wave velocity (v), porosity (n), wellbore radius (rw), storativity (S), or effective height of water column (He) in which each curve shows the changes of K caused by the increase of each input parameter while other parameter are unaltered. Figure 12. View largeDownload slide Model sensitivity to the estimated parameters: (a) curve smoothness of CSD caused by different length (L) of data, and the corresponding porosity (n) and storativity (S); (b) estimated values of hydraulic conductivity (K) by inputting different wave velocity (v), porosity (n), wellbore radius (rw), storativity (S), or effective height of water column (He) in which each curve shows the changes of K caused by the increase of each input parameter while other parameter are unaltered. From eqs (8) and (9), the amplification factors were determined by the wellbore radius rw, porosity n, wave velocity v, effective water volume height He (related to the thickness of the aquifer), storativity S and transmissivity T (or hydraulic conductivity K) if the other parameters were given, such as the bulk modulus of elasticity of water Ew, density of water ρ and gravitational acceleration g. As shown in Fig. 12(b), based on the data of the M7.8 earthquake on 2016 December 9 for Well X10, the estimated hydraulic conductivity using the model of Cooper et al. (1965) increases with an increase in rw, n and v and decreases with a decrease in S and He. Although the estimated hydraulic conductivity was affected by those input parameters, the magnitude of the estimated K did not change when the inputs were given as reasonable values. 4.3 Parameters estimated by different methods As shown in Fig. 8(a), the tidal response was clearer with the corrected water level than the raw water level, and a step-wise change induced by the M6.0 earthquake can also be identified. Many studies have employed the tidal response of the water level as a proxy for studying the hydraulic parameters of well-aquifer systems (Elkhoury et al.2006; Xue et al.2013; Shi & Wang 2016; Wang et al.2016), and the most commonly used method for this analysis is Hsieh's horizontal flow model (Hsieh et al.1987). When the aquifer is well confined, the permeability can be determined by the model using the negative phase shift between the water level and the tidal strain. However, Hsieh's model does not apply when the phase shift between the water level and the tidal strain is positive (Shi & Wang 2016). Roeloffs (1996) proposed that a positive phase shift results from the diffusion of pore pressure to the water table (the aquifer is not well confined), and the degree of confinement controls the frequency response of the pore pressure to tidal loading (Kinoshita et al.2015). Because hydraulic diffusivity is related to transmissivity and storativity, D = T/S, and the relationship between phase shift and hydraulic diffusivity can be written as follows:   $$p(z,\omega ) = B{K_u}\varepsilon (1 - {e^{ - (1 + i)z/\sqrt {2D/\omega } }})$$ (10)where p(z,ω) is the pore-pressure fluctuation at depth z, B is Skempton's coefficient, Ku is the bulk modulus of the saturated rock under undrained conditions, ɛ is the change in the volumetric strain and ω is the frequency of fluctuation. Roeloffs (1996) also mentioned that the response as a function of time (t) to a unit strain step (ɛ0) can be written as follows:   $$p(z,t) = - B{K_u}{\varepsilon _0}{\rm{erf}}(\sqrt {{z^2}/4Dt} ).$$ (11) By using the Venedikov's program of tidal analysis (Venedikov et al.2003), the phase shift of the water level in Well X10 was analysed. For the analysis, we set a data span of 28 d with a sliding step of one day. Then, we used the average value of the phase shifts of each month to determine the aquifer parameters. As shown in Fig. 13(a), the phase shift was obvious from October 2016 to January 2017. Based on Hsieh's horizontal flow model and Roeloffs’ vertical diffusion model, we can estimate the transmissivity and hydraulic diffusivity of aquifers (Shi & Wang 2016). To differentiate the parameters from the two models, we used subscript h to indicate the horizontal model and v for the vertical model. In addition, the water-level step-rise after M6.0 earthquake (marked by a dashed rectangle in Fig. 13a) was enlarged as shown in Fig. 13(b), and the hydraulic diffusivity was estimated according to eq. (11). Figure 13. View largeDownload slide Changes of water level and tidal phase shift in Xin 10 well during October 2016 to January 2017 (a). Water level and phase shift variations during October 2016 to January 2017 (b). Coseismic water level response following the 2016-12-08 M6.0 earthquake and fitting curve of the vertical model by using different hydraulic diffusivity values. Figure 13. View largeDownload slide Changes of water level and tidal phase shift in Xin 10 well during October 2016 to January 2017 (a). Water level and phase shift variations during October 2016 to January 2017 (b). Coseismic water level response following the 2016-12-08 M6.0 earthquake and fitting curve of the vertical model by using different hydraulic diffusivity values. The estimation results are listed in Table 2. Based on the coseismic response, the values of S and Th were estimated by the models mentioned in this study (Cooper et al.1965; Shih 2009); and based on the tidal response, the values of S and Th were estimated by Hsieh's horizontal flow model and Dh was calculated by Dh = Th/S. Then, based on Roeloffs’ vertical diffusion model, Dv was estimated by positive phase shifts and afterquake recovery and corresponding Tv was calculated by Tv = DvS. The results show that the value of storativity estimated from coseismic response was slightly smaller (less than an order of magnitude), although the value of transmissivity was much higher (approximately two orders of magnitude) than the value inferred from tidal response. The value of vertical diffusivity Dv estimated from tidal response was similar to the value from the afterquake recovery. Table 2. Aquifer parameters estimated from different data.     Transmissivity  Hydrulic diffusivity  Responses  Storativity (10−6)  Th (m2 d−1)  Tv (m2 d−1)  Dh (m2 d−1)  Dv (m2 d−1)  Coseismic response  10–20  320–800  —  2–5 × 107  —  Tidal response  70–120a  2–14a  0.002–0.004b  0.2–2 × 105a  30.83b  Stain step response  —  —  —  —  20–50      Transmissivity  Hydrulic diffusivity  Responses  Storativity (10−6)  Th (m2 d−1)  Tv (m2 d−1)  Dh (m2 d−1)  Dv (m2 d−1)  Coseismic response  10–20  320–800  —  2–5 × 107  —  Tidal response  70–120a  2–14a  0.002–0.004b  0.2–2 × 105a  30.83b  Stain step response  —  —  —  —  20–50  aFrom the negative phase shift in Fig. 13(a). bFrom the positive phase shift in Fig. 13(b). Notes: The em-dash sign (—) means it is not applicable. View Large The higher estimated value in transmissivity by the seismic wave method is not caused by an inappropriate model or error in estimation but rather to the different frequency response of the well-aquifer system. As with the barometric response (Fig. 8b), the response coefficient varies with the frequency of loading. For Well X10 aquifer system, it is likely that the higher the loading frequency, the greater the transmissivity (or hydraulic conductivity). Under periodic loading, especially in the short term, such as a seismic wave with a period less than one minute, the water in the aquifer might not flow far. As in Darcy's law (K = v/I), for a given hydraulic gradient I, periodic loading in the short term will speed up the flow of water (v), which will increase the effective hydraulic conductivity, K. Therefore, the transmissivity estimated from the coseismic response will be higher than the value from the tidal response. In addition to the response from coseismic and barometric loading, this mechanism is also reflected in the response to wind. From the water level (raw or corrected) shown in Fig. 8(a), fluctuations caused by short-term strong winds are also obvious. In addition, the parameter values estimated from the horizontal flow model were higher than the values from the vertical flow model (Table 2), which may indicate that the water flow in Well X10 aquifer is mainly horizontal. 4.4 Scale effect when estimating hydraulic parameters Accurate observations are important for obtaining reliable results. The Rayleigh wave timescale is approximately tens of seconds. Therefore, the sampling rate of the seconds scale is critical, followed by the resolution on the millimetre scale. As shown in Table 1, the water-level amplitudes induced by the seismic wave range from 3.5 to 30 mm, which requires the resolution of the pressure sensor be at least 1 mm. Moreover, the precision and accuracy of the sensor must be higher. In this study, the pressure sensors have a resolution of 1 mm and sampling rate of 1 s, which provide a guarantee for the complete record of water-level fluctuation. The sensor is located approximately 1 m below the surface of the water, with a precision of 0.02 per cent and accuracy of 0.05 per cent. The F.S. measurement error (repeatability of successive measurements under the same conditions) is approximately 0.2 mm, and the accuracy (the difference between the measured value and the true) is 0.5 mm. Because we are concerned with the relative changes in water level, the measurement error of 0.2 mm meets our research needs. To obtain accurate estimates of hydraulic conductivity, water must flow within the well-aquifer system. The water-level fluctuation in the well was caused by the periodic loading of seismic waves, and the motion of the water column created a periodic discharge from the aquifer into the well accompanied by drawdown, which is approximated by Cooper et al. (1965)   $$s = \frac{{\omega r_w^2{x_0}}}{{2T}}[Ker(\alpha )\cos (\omega t) - Kei(\alpha )\sin (\omega t)]$$ (9)where $$\alpha = r\sqrt {\omega S/T}$$ (dimensionless); x0 is the amplitude of the water-level fluctuation; r is the radial distance from the centre of well; and $$\omega r_w^2{x_0}/2T$$ represents the effect of the hydraulic gradient, which mainly affects the amplitude of drawdown. Ker(α)cos (ωt) − Kei(α)sin (ωt) represents the effective water storage of the aquifer; thus, it determined the range that produces significant drawdown. If we set S = 10−5, Tτ = 0.1 m2 and rw = 6 cm, the region of significant drawdown is approximated by a radius of 100 m. Hence, the parameters we estimate from coseismic water-level fluctuation can represent the hydraulic conductivity/transmissivity at a scale of hundreds of metres. In addition to the amplitude of water-level response to surface motion, the phase angle is as follows:   $$\eta = \arctan \left[\frac{{\omega r_w^2gKer({\alpha _w})}}{{2T{H_e}({\omega ^2} - \omega _w^2)}}\right]$$ (10)where $${\omega _w} \approx \sqrt {g/{H_e}}$$, which is the resonant frequency of a well. According to eq. (10), the phase ranges from −90° to 90° and decreases with an increase in frequency. The phase is negative when the frequency is less than the natural frequency and positive when the natural frequency is greater than the natural frequency; in addition, the phase is reversed when the frequency is equal to the natural frequency. As in the barometric efficiency/response method, the seismic wave and tide response methods that estimate hydraulic conductivity are also frequency dependent. Additionally, the hydraulic conductivity may vary over an order of magnitude on a small scale because it is dependent on the scale at which it is measured. Thus, discussing the scale effect on the method is necessary. Because the Cooper et al. (1965) and Hsieh et al. (1987) models were derived from the model of Theis (1935), their assumptions and solutions were similar. Here, if we assume that the amplitudes of periodic water-level fluctuation are the same, for example, 10 cm with a well radius of 6 cm, the amplitude of the discharge induced by the seismic wave will be on a dozens of cubic metres scale, whereas only dozens of litres for the tidal response would induce discharge. The former is more similar to the discharge produced in the field pumping tests. Because of the frequency scale, the model of Hsieh et al. is more suitable for determining aquifer transmissivities less than 1 m2 d−1; otherwise, the estimated value will be smaller than the actual value. The model of Cooper et al. is more suitable for aquifers with larger transmissivity. Hence, the transmissivity results, Th, inferred from the tidal response are smaller than the actual values as shown in Table 2. In addition, the large hydraulic diffusivity values of Dh may not represent the actual value of the whole aquifer system because of (1) the several hundred metres scale and periodic loading, and (2) the assumption that the release of groundwater from storage is instantaneous because of the drop in water head. 5 CONCLUSIONS This paper presents a method for estimating the hydraulic parameters of confined aquifers using water-level oscillations induced mainly by Rayleigh waves. The methodology is demonstrated using the coseismic responses of water levels in Well X10 to distant earthquakes and correlations with seismic waves that were analysed using spectral methods. The water-level oscillations in Well X10 are correlated with seismic waves in both their amplitude and frequency, and the aquifer hydraulic properties inferred from the seismic wave response of Well X10 yield a hydraulic conductivity of K = 1 × 10−3 m s−1 and a storativity of S = 15 × 10−6. The method proposed in this study is especially useful for wells that must be observed continuously and for which a clear tidal signal cannot be recorded. Acknowledgements The original water level and seismic wave data supporting the work presented in this paper are from the monitoring centre of the Earthquake Agency of the Xinjiang Uygur Autonomous Region. These data were extremely important to the study and the knowledge gained. We thank Todd C. Rasmussen and an anonymous reviewer for their constructive reviews and the help provided by the Editor. This research was supported by the National Natural Science Foundation of China (41502239 to Sun XL, U1602233 and 41602266 to Shi ZM), the National Key R&D Program of China (2017YFC1500502 to Sun XL) and the Research Grant from the Institute of Crustal Dynamics, China Earthquake Administration (ZDJ2017-27 to Sun XL). REFERENCES Allègre V., Brodsky E.E., Xue L., Nale S.M., Parker B.L., Cherry J.A., 2016. Using earth-tide induced water pressure changes to measure in situ permeability: a comparison with long-term pumping tests, Water Resour. Res. , 52( 4), 3113– 3126. https://doi.org/10.1002/2015WR017346 Google Scholar CrossRef Search ADS   Bendat J.S., Piersol A.G., 2000. Random Data, Analysis and Measurement Procedures , John Wiley and Sons. Brodsky E.E., Roeloffs E., Woodcock D., Gall I., Manga M., 2003. A mechanism for sustained groundwater pressure changes induced by distant earthquakes, J. geophys. Res. , 108( B8), doi:10.1029/2002JB002321. https://doi.org/10.1029/2002JB002321 Candela T., Brodsky E.E, Marone C., Elsworth D., 2014. Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing, Earth plant. Sci. Lett. , 392, 279– 291. https://doi.org/10.1016/j.epsl.2014.02.025 Google Scholar CrossRef Search ADS   Chapuis R.P., 1992. Similarity of internal stability criteria for granular soils, Can. Geotech. J. , 29( 4), 711– 713. https://doi.org/10.1139/t92-078 Google Scholar CrossRef Search ADS   Clark W.E., 1967. Computing the barometric efficiency of a well, J. Hydraul. Div. , 93( 4), 93– 98. Cooper H.H., Bredehoeft J.D., Papadopulos I.S., Bennett R.R., 1965. The response of well-aquifer systems to seismic waves, J. geophys. Res. , 70( 16), 3915– 3926. https://doi.org/10.1029/JZ070i016p03915 Google Scholar CrossRef Search ADS   Cooper H.H., Jacob C.E., 1946. A generalized graphical method for evaluating formation constants and summarizing well-field history, EOS, Trans. Am. geophys. Un. , 27( 4), 526– 534. https://doi.org/10.1029/TR027i004p00526 Google Scholar CrossRef Search ADS   Dawson K.J., Istok J.D., 1991. Aquifer Testing: Design and Analysis of Pumping and Slug Tests , CRC Press. Duhamel P., Vetterli M., 1990. Fast Fourier transforms: a tutorial review and a state of the art, Signal Process. , 19( 4), 259– 299. https://doi.org/10.1016/0165-1684(90)90158-U Google Scholar CrossRef Search ADS   Elkhoury J.E., Brodsky E.E., Agnew D.C., 2006. Seismic waves increase permeability, Nature , 441( 7097), 1135– 1138. https://doi.org/10.1038/nature04798 Google Scholar CrossRef Search ADS PubMed  Elkhoury J.E., Niemeijer A., Brodsky E.E., Marone C., 2011. Laboratory observations of permeability enhancement by fluid pressure oscillation of in situ fractured rock, J. geophys. Res. , 116( B2), doi:10.1029/2010JB007759. https://doi.org/10.1029/2010JB007759 Freeze R.A., Cherry J.A., 1979. Groundwater , Prentice-Hall. Hsieh P.A., Bredehoeft J.D., Farr J.M., 1987. Determination of aquifer transmissivity from earth tide analysis, Water Resour. Res. , 23( 10), 1824– 1832. https://doi.org/10.1029/WR023i010p01824 Google Scholar CrossRef Search ADS   Jacob C.E., 1940. On the flow of water in an elastic artesian aquifer, EOS, Trans. Am. geophys. Un. , 21( 2), 574– 586. https://doi.org/10.1029/TR021i002p00574 Google Scholar CrossRef Search ADS   Jiao J.J., Rushton K.R., 1995. Sensitivity of drawdown to parameters and its influence on parameter estimation for pumping tests in large-diameter wells, Ground Water , 33( 5), 794– 800. https://doi.org/10.1111/j.1745-6584.1995.tb00026.x Google Scholar CrossRef Search ADS   Kinoshita C., Kano Y., Ito H., 2015. Shallow crustal permeability enhancement in central Japan due to the 2011 Tohoku earthquake, Geophys. Res. Lett. , 42( 3), 773– 780. https://doi.org/10.1002/2014GL062792 Google Scholar CrossRef Search ADS   Ma L., Xu Y.S., Shen S.L., Sun W.J., 2014. Evaluation of the hydraulic conductivity of aquifers with piles, Hydrogeol. J. , 22( 2), 371– 382. https://doi.org/10.1007/s10040-013-1068-y Google Scholar CrossRef Search ADS   Manga M., Beresnev I., Brodsky E.E., Elkhoury J.E., Elsworth D., Ingebritsen S.E., Wang C.Y., 2012. Changes in permeability caused by transient stresses: field observations, experiments, and mechanisms, Rev. Geophys. , 50( 2), doi:10.1029/2011RG000382. https://doi.org/10.1029/2011RG000382 Rasmussen T.C., Crawford L.A. 1997. Identifying and removing barometric pressure effects in confined and unconfined aquifers, Groundwater , 35( 3), 502– 511. https://doi.org/10.1111/j.1745-6584.1997.tb00111.x Google Scholar CrossRef Search ADS   Roeloffs E., 1996. Poroelastic techniques in the study of earthquake-related hydrologic phenomena, Adv. Geophys. , 37, 135– 195. Google Scholar CrossRef Search ADS   Sethi R., 2011. A dual-well step drawdown method for the estimation of linear and non-linear flow parameters and wellbore skin factor in confined aquifer systems, J. Hydrol. , 400( 1), 187– 194. https://doi.org/10.1016/j.jhydrol.2011.01.043 Google Scholar CrossRef Search ADS   Shi Z., Wang G., 2016. Aquifers switched from confined to semiconfined by earthquakes, Geophys. Res. Lett. , 43( 21), 11 166–11 172. Google Scholar CrossRef Search ADS   Shi Z., Wang G., Manga M., Wang C., 2015. Mechanism of co-seismic water level change following four great earthquakes—insights from co-seismic responses throughout the Chinese mainland, Earth planet. Sci. Lett. , 430, 66– 74. https://doi.org/10.1016/j.epsl.2015.08.012 Google Scholar CrossRef Search ADS   Shih D.C.F., 2009. Storage in confined aquifer: spectral analysis of groundwater responses to seismic Rayleigh waves, J. Hydrol. , 374( 1), 83– 91. https://doi.org/10.1016/j.jhydrol.2009.06.002 Google Scholar CrossRef Search ADS   Shih D.C.F., Wu Y.M., Chang C.H., 2013. Significant coherence for groundwater and Rayleigh waves: evidence in spectral response of groundwater level in Taiwan using 2011 Tohoku earthquake, Japan, J. Hydrol. , 486, 57– 70. https://doi.org/10.1016/j.jhydrol.2013.01.013 Google Scholar CrossRef Search ADS   Stockwell R.G., Mansinha L., Lowe R.P., 1996. Localization of the complex spectrum: the S transform, IEEE Trans Signal Process , 44( 4), 998– 1001. Google Scholar CrossRef Search ADS   Sun X, Wang G, Yang X., 2015. Coseismic response of water level in Changping well, China, to the Mw 9.0 Tohoku earthquake, J. Hydrol. , 531, 1028– 1039. https://doi.org/10.1016/j.jhydrol.2015.11.005 Google Scholar CrossRef Search ADS   Theis C.V., 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, EOS, Trans. Am. geophys. Un. , 16( 2), 519– 524. https://doi.org/10.1029/TR016i002p00519 Google Scholar CrossRef Search ADS   Venedikov A.P., Arnoso J., Vieira R., 2003. VAV: a program for tidal data processing, Comput. Geosci. , 29 ( 4), 487– 502. https://doi.org/10.1016/S0098-3004(03)00019-0 Google Scholar CrossRef Search ADS   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( 8), 5923– 5937. https://doi.org/10.1002/2016WR018893 Google Scholar CrossRef Search ADS   Wang C.-Y., 2007. Liquefaction beyond the Near Field, Seismol. Res. Lett. , 78( 5), 512– 517. https://doi.org/10.1785/gssrl.78.5.512 Google Scholar CrossRef Search ADS   Wang C.-Y., Manga M., 2010. Hydrologic responses to earthquakes and a general metric, Geofluids. , 10( 1–2), 210– 216. Xue L., Li H.B., Brodsky E.E., Xu Z.Q., Kano Y., Wang H., Yang G., 2013. Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone, Science , 340( 6140), 1555– 1559. https://doi.org/10.1126/science.1237237 Google Scholar CrossRef Search ADS PubMed  © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Estimating the hydraulic parameters of a confined aquifer based on the response of groundwater levels to seismic Rayleigh waves

, Volume 213 (2) – May 1, 2018
12 pages

/lp/ou_press/estimating-the-hydraulic-parameters-of-a-confined-aquifer-based-on-the-KhLiDqkk1o
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy036
Publisher site
See Article on Publisher Site

### Abstract

Summary Groundwater flow models implemented to manage regional water resources require aquifer hydraulic parameters. Traditional methods for obtaining these parameters include laboratory experiments, field tests and model inversions, and each are potentially hindered by their unique limitations. Here, we propose a methodology for estimating hydraulic conductivity and storage coefficients using the spectral characteristics of the coseismic groundwater-level oscillations and seismic Rayleigh waves. The results from Well X10 are consistent with the variations and spectral characteristics of the water-level oscillations and seismic waves and present an estimated hydraulic conductivity of approximately 1 × 10−3 m s−1 and storativity of 15 × 10−6. The proposed methodology for estimating hydraulic parameters in confined aquifers is a practical and novel approach for groundwater management and seismic precursor anomaly analyses. Permeability and Porosity Hydrogeophysics, Hydrology, Fourier analysis, Time series analysis 1 INTRODUCTION Understanding hydrological properties is crucial for many applications, such as groundwater resource management, subsurface solute transport assessments, hydraulic fracturing and mining. Transmissivity (T) and storativity (S) are the two key parameters that must be accurately determined in many of these applications. The methods for estimating T and S include laboratory experiments, field tests and water-level responses in wells caused by natural perturbations. Laboratory methods analyse the physical properties (such as the porosity, grain size and compression coefficient) of the rock or soil samples (Chapuis 1992; Ma et al.2014). Field tests are conducted by pumping a well and then observing water-level responses in the pumped well or in nearby observation wells. Data from these tests are interpreted using a range of analytic solutions, including the Theis, corrected Theis and Hantush models (Theis 1935; Jacob 1940; Cooper & Jacob 1946). The complex geometry, boundary conditions and heterogeneity parameterizations have been incorporated into numerous analytical pump test solutions (Dawson & Istok 1991; Jiao & Rushton 1995; Sethi 2011). Water-level responses to natural perturbations include barometric-pressure fluctuations, surface loading and earth tides (Hsieh et al.1987; Allègre et al.2016). Of these three methods, laboratory methods may provide inaccurate parameter values because of the natural heterogeneity and artefacts associated with the specific testing protocol. Although field testing is common, it can be expensive, and continuous observations of monitoring wells may not be feasible. The response to natural perturbations is becoming more widely used because real-time monitoring systems are increasingly adopted and appropriate methods for determining hydraulic parameters have been developed (Xue et al.2013; Wang et al. 2016; Shi & Wang 2016), especially the responses of well-water levels to earth tides to determine the aquifer hydraulic conductivity. However, this method is hindered because it requires a clear tidal signal in the observation wells. Water-level responses to earthquakes represent another method for estimating aquifer properties (Shih 2009; Sun et al.2015; Wang et al.2016). The main coseismic groundwater-level response types are step change and oscillation change (Shi et al.2015). Coseismic step changes have been used to estimate permeability (Wang et al.2016), and oscillation changes have been used to estimate specific storage (Shih 2009; Sun et al.2015). However, to the best of our knowledge, a method of estimating aquifer hydraulic properties from coseismic responses has not been proposed. Therefore, we present a method that compares the correlation between the spectral characteristics of coseismic water-level oscillation and seismic Rayleigh waves, which we then use to identify the hydraulic properties of a well in a confined aquifer. 2 BACKGROUND AND OBSERVATIONS Well X10 is located to the south of Urumqi City in the Xinjiang Autonomous Region in western China at 43.70° N and 87.62° E and an elevation of 1056 m amsl. The well hole was drilled in September 1980. Well X10 is at the intersection of the Liushugou-Hongyanchi thrust fault and its derived fault (Fig. 1). Groundwater discharges through springs at this location provide stable year-round yields of 20–30 l s−1. Figure 1. View largeDownload slide Geological setting, location of Well X10 (blue triangle) and epicentres of the 10 earthquakes (beach balls) that occurred between November and December 2016 Figure 1. View largeDownload slide Geological setting, location of Well X10 (blue triangle) and epicentres of the 10 earthquakes (beach balls) that occurred between November and December 2016 As shown in Fig. 2, the total depth of Well X10 is 28 m, the diameter is 146 mm at depths to 9.2 m, the diameter is 130 mm at depths from 9.2 to 13.04 m, and the diameter is 110 mm at depths from 13.04 m to the bottom. Core samples collected between 16 and 24 m showed extensive fracturing; therefore, a screen filter was installed in this zone. The aquifer in the Permian Hongyanchi group is 8 m thick, consists of sandstone and gravel and has a static water level of approximately 1 m below ground surface. Figure 2. View largeDownload slide Aquifer lithology and wellbore structure of Well X10. Figure 2. View largeDownload slide Aquifer lithology and wellbore structure of Well X10. Starting in October 2016, the water level was monitored using an SWY-II type digital water-level metre, which was developed and produced by the Institute of Crustal Dynamics, China Earthquake Administration. The water-level sensor has a range of 0–50 m, resolution of 1 mm, sampling interval of 1 s, precision of 0.02 per cent and accuracy of 0.05 per cent F.S. Ten earthquakes during November and December 2016 caused the water level to fluctuate in Well X10 (Fig. 1). The water-level changes in a well induced by far-field earthquakes were strongly correlated to the dynamic strains. Wang (2007) introduced a physical quantity, the seismic energy density e (J m−3), which is a function of the earthquake magnitude (M) and epicentral distance, r (km): log r = 0.48M − 0.33 log e – 1.4. A comparison with other well sites and other earthquake-related hydrological phenomena, we plotted the earthquake events that induced water-level changes using the magnitude versus epicentral distance to Well X10 (Fig. 3). A global data set of water-level responses compiled by Wang & Manga (2010) were also plotted for comparison. The seismic energy densities were significantly lower in Well X10 than in the global data set, which were mainly within 10−4–10−7 J m−3, indicating that Well X10 is highly sensitive to seismic waves. Figure 3. View largeDownload slide Seismic energy density as a function of epicentral distance and earthquake magnitude. Lines of constant seismic energy density are given in J m−3, and the red squares and grey circles correspond to Well X10 and a global data set compiled by Wang & Manga (2010), respectively. Figure 3. View largeDownload slide Seismic energy density as a function of epicentral distance and earthquake magnitude. Lines of constant seismic energy density are given in J m−3, and the red squares and grey circles correspond to Well X10 and a global data set compiled by Wang & Manga (2010), respectively. The fluctuations in Well X10 are compared with vertical displacement waveforms observed by the SLJ-100 seismic recoding sensor at station SMG, which is located 10 km northeast of Well X10. Fig. 4 shows a comparison of the coseismic responses of the water levels in Well X10 with the vertical displacement waveforms at station SMG caused by the M 7.8 Solomon Islands earthquake (10.665° S, 161.335° E), which occurred at 01:38 China Standard Time (CST) on 2016 December 9. Figure 4. View largeDownload slide Water level changes in the Well X10 (a) and vertical displacement at SMG (b) due to the 2016 December 9 M7.8 Solomon Islands earthquake and time–frequency images of water level (c) and displacement (d). The colour bar shows the amplitude of water level or displacement. Figure 4. View largeDownload slide Water level changes in the Well X10 (a) and vertical displacement at SMG (b) due to the 2016 December 9 M7.8 Solomon Islands earthquake and time–frequency images of water level (c) and displacement (d). The colour bar shows the amplitude of water level or displacement. To better illustrate the similarities between the two waveforms, the S-transform method (Stockwell et al. 1996) was used to calculate the time–frequency characteristics as presented in Fig. 4(c) (water-level oscillation in Well X10) and Fig. 4(d) (seismic waves at the SMG station). Time is on the horizontal axis, frequency is on the vertical axis, and the colours represent response amplification. Note that the similarities in the shapes and frequencies of the water-level oscillations in Well X10 and the vertical component of the seismic waves at the SMG station were caused by the M7.8 Solomon Islands earthquake. The S-waves and surface waves were also clearly detected. A comparison of the time–frequency characteristics in Figs 4(c) and (d) show that the water and seismic waves also have a wide frequency range of 0.03–0.06 Hz (15–35 s) and last approximately 20 min. For the surface wave, the low frequency and long period signals were first observed and followed by high frequency and short period signals. The comparison of the amplifications shows that the water and seismic waves reached a range of values at a period of approximately 20–40 s, which was the main period of the surface wave (Rayleigh wave). 3 ESTIMATION OF HYDRAULIC PARAMETERS The analysis presented above highlights the similarity between the shapes and frequencies of the water-level oscillations and the vertical component seismic waves, especially for the surface wave period. Cooper et al. (1965) investigated the responses of aquifer systems to vertical surface movements and stated that the water-level responses were induced by both the vertical surface movements and pore-pressure fluctuations in the aquifer (Fig. 5), and their results also demonstrated that the response amplifications varied with the signals at different periods and were closely correlated with the hydraulic parameters (e.g. transmissivity and storativity). Figure 5. View largeDownload slide Idealized representation of an open well in which water level fluctuations are caused by the oscillation of artesian pressure and vertical motion of land surface. Where p = ρg(H + z) + p0 sin (ωt − η) represents the pressure fluctuation in aquifer with a frequency ω and phase lag η; h = h0 sin (ωt) means the water-level fluctuation caused by pressure fluctuation; and s = s0 sin (ωt − ζ) represents the fluctuation caused by surface land motion with a phase lag ζ. Figure 5. View largeDownload slide Idealized representation of an open well in which water level fluctuations are caused by the oscillation of artesian pressure and vertical motion of land surface. Where p = ρg(H + z) + p0 sin (ωt − η) represents the pressure fluctuation in aquifer with a frequency ω and phase lag η; h = h0 sin (ωt) means the water-level fluctuation caused by pressure fluctuation; and s = s0 sin (ωt − ζ) represents the fluctuation caused by surface land motion with a phase lag ζ. As a result, hydraulic parameters can be estimated by analysing the water levels during the surface wave period and the magnitudes of the ground motions (Sun et al.2015). Shih (2009) used the spectral analysis method to investigate the water response to seismic surface waves and estimated the storativity from spectral correlations of the water levels and the vertical displacement of the earth's surface. These methods are used to estimate aquifer hydraulic parameters by analysing the characteristics of the water-level oscillations in Well X10 in response to 10 earthquakes from around the world. 3.1 Storativity and specific storage The storativity (S, dimensionless) is defined as S = d·Ss, where d is the thickness of aquifer and Ss is defined as the volume of water released from a unit volume of porous aquifer when there is a unit decline in the head (Freeze & Cherry 1979), in the dimensions of length−1. The amplitude of the response of water levels to ground motions is closely correlated with aquifer thickness, storativity and surface wave number (Shih 2009):   $$h(t) = \frac{{0.918d{k_x}}}{S}w(t),$$ (1)where h(t) is the water level (dimensions of length), w(t) is the ground displacement (dimensions of length), and kx is the wave number on the surface. The amplitudes of the responses to signals at different periods were calculated using Fourier transforms (Duhamel & Vetterli 1990):   $$\begin{array}{@{}l@{}} {Z_h}(\omega ) = \int_{{ - \infty }}^{\infty }{{h(t){{\rm{e}}^{ - 2\pi i\omega }}^t{\rm{d}}t}}\\ {Z_w}(\omega ) = \int_{{ - \infty }}^{\infty }{{w(t){{\rm{e}}^{ - 2\pi i\omega }}^t{\rm{d}}t}}, \end{array}$$ (2)where Zh(ω) and Zw(ω) are the Fourier transforms of h(t) and w(t), respectively, ω is the frequency (dimensions of angular per time), t is time, and i = $$\sqrt { - 1} \$$. The time-series and frequency spectra of the water level and vertical displacement induced by the M7.8 Solomon Islands earthquake are plotted in Fig. 6. Compared with the waves at other periods, the response amplitude at the 20–30 s period was relatively large. Figure 6. View largeDownload slide Time-series and its corresponding spectral of water level and vertical displacement at the surface wave period. Figure 6. View largeDownload slide Time-series and its corresponding spectral of water level and vertical displacement at the surface wave period. The auto-spectral density (PSD) is used to detect the strong signal in time-series, whereas the cross-spectral density (CSD) and coherence measure the intensity of specific target signals between two time-series. The 95 per cent confidence interval of the PSD and non-zero coherence level for the CSD are also evaluated to identify significant components in the frequency domain (Shih et al.2013). Aquifer storativity was calculated from the spectral correlations of the water level and vertical displacement as follows (Shih 2009):   $$S = \frac{{1.836\pi d}}{\lambda }{\left( {\frac{{{S_{hh}}}}{{{S_{ww}}}}} \right)^{ - 1/2}}$$ (3)  $$S = \frac{{1.836\pi d}}{\lambda }{\left( {\frac{{{S_{hh}}}}{{{S_{hw}}}}} \right)^{ - 1}}$$ (4)  $$S = \frac{{1.836\pi d}}{\lambda }{\left( {\frac{{{S_{wh}}}}{{{S_{ww}}}}} \right)^{-1}}$$ (5)where λ is the wave length (dimensions of length), which can be determined from the seismic velocity and period; Shh and Sww are the PSD of the water level h(t) and the vertical displacement of ground motion w(t), respectively; and Swh and Shw are the CSDs of h(t) and w(t), respectively. PSD and CSD can be found using a finite length (L) of the observation data (Bendat & Piersol 2000):   \begin{eqnarray} {S_{hh}}(\omega) &=& \mathop {\lim }\limits_{L \to \infty } \frac{2}{L}E\left[ {|{Z_h}(\omega ,L){|^2}} \right]\nonumber\\ {S_{ww}}(\omega) &=& \mathop {\lim }\limits_{L \to \infty } \frac{2}{L}E\left[ {|{Z_w}(\omega ,L){|^2}} \right]\\ {S_{hw}}(\omega) &=& \mathop {\lim }\limits_{L \to \infty } \frac{2}{L}E\left[ {{Z_w}(\omega ,L) \cdot {Z_w}(\omega ,{L_r})} \right],\nonumber \end{eqnarray} (6)where the expected value E is the average value in a sampling period. The PSD and CSD of the water level and vertical displacement induced by the M7.8 Solomon Islands earthquake are plotted in Figs 7(a) and (b), respectively. The PSD with the 95 per cent confidence interval has lower and the upper extremes of 0.7364 and 3.8797 cm2 Hz−1, respectively (Fig. 7a). The non-zero coherence significance level is 0.68 (Fig. 7b). The plots indicate that the spectral densities of the water level and vertical displacement were largest and highly correlated in the 0.03–0.04 frequency range and could be used to estimate the storativity (Shih et al.2013). Figure 7. View largeDownload slide Autospectral density of water level and vertical displacement (a), and cross-spectral density (b), the short dash lines indicate the 95 per cent non-zero coherence level. Figure 7. View largeDownload slide Autospectral density of water level and vertical displacement (a), and cross-spectral density (b), the short dash lines indicate the 95 per cent non-zero coherence level. The water bearing layer is at depths of 16–24 m in Well X10, and its thickness is d = 8 m. The wave velocity v in the formula for the wavelength calculation λ = v/ω is obtained by dividing the distance to the earthquake epicentre by the difference between the arrival time of surface wave and the start time of water-level oscillation. The frequency ω was the highest point of the PSD in Fig. 7(a). According to eqs (3)—(6) and the spectral correlations of the water level and vertical displacement induced by the M7.8 Solomon Islands earthquake, the average estimated storativity in the Well X10 was 19.2 × 10−6. 3.2 Aquifer porosity Porosity (n) is the pore volume per unit aquifer volume. Jacob (1940) derived an equation that relates specific storage to porosity:   $$n = \frac{{{B_e}{S_s}}}{{\beta \rho }},$$ (7) where Be is the barometric efficiency (dimensionless), β = 4.6 × 10−10 m2 N−1 is the water compressibility and ρ = 103 kg m−3 is the fluid density. The barometric efficiency (Be) was estimated using Clark's method (Clark 1967) and the regression deconvolution method (Rasmussen & Crawford 1997). As shown in Fig. 8(a), the water level in Well X10 was influenced by the barometric pressure. We calculated the barometric efficiency based on the changes (differential values) in water level and barometric pressure and obtained the aquifer response from both the differential values and detrended values. The results show that the barometric efficiency and the aquifer response decrease gradually as the time interval increases when using the differential values, and the aquifer response is stable when using the detrended values (Fig. 8b). The water-level response to barometric pressure based on the differential values could reveal the short-term response, and the detrended values could reveal a long-term response. Therefore, we concluded that the barometric efficiency of Well X10 varies with the barometric-pressure period. Figure 8. View largeDownload slide Comparison of water level and barometric pressure (a) and the barometric efficiency/response corresponding to different time (b) in Well X10 between 2016 November 1 and December 31. Be is the barometric efficiency with different time interval, and Br is the barometric response with different time lag. The segment of Br corresponding to the time between 1 and 500 min was calculated by minute values, and the segment of Br corresponding to the time between 60 and 10 000 min was calculated by hourly values. Figure 8. View largeDownload slide Comparison of water level and barometric pressure (a) and the barometric efficiency/response corresponding to different time (b) in Well X10 between 2016 November 1 and December 31. Be is the barometric efficiency with different time interval, and Br is the barometric response with different time lag. The segment of Br corresponding to the time between 1 and 500 min was calculated by minute values, and the segment of Br corresponding to the time between 60 and 10 000 min was calculated by hourly values. The time intervals used in the method were selected such that the barometric-pressure-independent water-level change was as small as possible compared with the barometric-pressure change. Focusing the time-interval selection around the correlation coefficient maximizes the value of the barometric-pressure change while minimizing the barometric-pressure-independent water-level change. As shown in Fig. 8(b), the barometric efficiency value was 0.5036, which was the mean of the barometric efficiencies with a correlation coefficient larger than 0.75. In combination with the specific storage Ss, the porosity of the aquifer can be estimated. 3.3 Transmissivity and hydraulic conductivity The response amplitude of the water-level oscillations induced by seismic waves is closely related to the aquifer parameters, borehole structure and seismic wave periods. In addition to causing vertical fluctuations in the surface medium, seismic waves can also cause pore-pressure fluctuations in the aquifer. As shown in Fig. 5, the oscillations of the well-water level are affected by both the vertical motions of the earth's surface and the pore-pressure changes, which are in turn amplified by the water levels. Therefore, the response amplitude of the water-level oscillations is affected by both the pore-pressure fluctuations in well-aquifer systems and vertical motions of the earth's surface. In a confined well-aquifer system, the response of well-water levels to the pore pressure and vertical surface displacement can be represented by an amplification factor (or amplitude ratio) (Cooper et al.1965):   \begin{eqnarray} \left. \begin{array}{@{}l@{}} A = {\left[{\left(1 - \displaystyle\frac{{\pi r_w^2}}{{T\tau }}Kei{\alpha _w} - \frac{{4{\pi ^2}{H_e}}}{{{\tau ^2}g}} \right)^2} + {\left(\frac{{\pi r_w^2}}{{T\tau }}Ker{\alpha _w}\right)^2}\right]^{ - 1/2}}\\ A' = A \cdot \displaystyle\frac{{4{\pi ^2}{H_e}}}{{{\tau ^2}g}},{\alpha _w} = {r_w}{(\omega S/T)^{1/2}} \end{array} \right\}\nonumber\\ \end{eqnarray} (8)where A (dimensionless) is the water-level amplification factor of the aquifer's pore-pressure fluctuations, A’ (dimensionless) is the water-level amplification factor of the vertical surface displacements, rw is the wellbore radius(dimensions of length), He is the effective height of water column (dimensions of length), T is the transmissivity (dimensions of square length per time), S is the storativity, τ is the period of seismic wave (dimensions of time), ω is the angular frequency, g is the gravitational acceleration, and Kei et Ker are the imaginary and real parts of the zero-order Kelvin functions. If the height of the water-level oscillations induced by fluctuations of the aquifer's pore pressure is h and the height of the water-level oscillations induced by the vertical surface displacement is s, then the total height of the water-level oscillations is x = h + s. If R is the ratio of the height induced by pore pressure to the height induced by the vertical displacement and M is the amplitude ratio of the total height of the water-level to the vertical displacement, then   $$R = \frac{h}{s} = \frac{{2.7{E_w}}}{{\rho gnv\tau }},M = \frac{x}{s} = \left(\frac{x}{h}\right) \cdot \left(\frac{h}{s}\right) = A \cdot R.$$ (9) Water at 18° centigrade has a bulk modulus of elasticity of Ew = 2.2 × 109 N m−2 and a density ρ of 103 kg m−3 (Cooper et al.1965); the gravitational acceleration is g = 9.81 m s−2; v and τ are the wave velocity and wave period, respectively; and n is the aquifer porosity. R and M were calculated using eqs (8) and (9) and the observed amplitudes of the water level and vertical displacement. If the ratio of the amplification factors A/A’ = R’, the ratio of the fluctuations R has an inverse relationship with the seismic wave period τ, although R’ has a positive quadratic relationship with τ. The calculation of R requires the wave velocity v, which can be obtained from the distance from the earthquake epicentre divided by the time difference between the arrival time of surface wave and the start time of the water-level oscillations. The period τ can be determined from the time–frequency diagram in Fig. 6. The effective water column height He required by the R’ calculation is He = H + 3d/8 (Cooper et al.1965), where H is the well's water column height above the upper confining bed. Because the upper confining bed in Well X10 is at a depth of 16 m and the depth to the water table is 1 m, the thickness of the well's water column is 15 m; therefore, He = H + 3d/8≈18 m. The M values at every period τ can be obtained statistically from the time–frequency values, and the respective values of A and A’ can be calculated with their R values (Fig. 9). The values represented by green circles in Fig. 9 were calculated using the spectral amplitude ratio of water level and displacement shown in Fig. 6. The magenta solid lines in Fig. 9 were the theoretical curves that fit the spectral amplitude ratio. To identify the optimal hydraulic conductivity (K) that can best fit the observed data, one thousand groups of theoretical values (A and A’) determined by eq. (8) were calculated with different values of K (from 0.001 to 1 cm s−1), and the standard deviations of each group were obtained. Then, the optimal K was selected from the group with the smallest standard deviation. Figure 9. View largeDownload slide Observed (circles) and estimated (black solid line) amplification factors based on the water-level oscillations in the Well X10: (a) due to pore pressure fluctuations; (b) due to vertical displacement. Figure 9. View largeDownload slide Observed (circles) and estimated (black solid line) amplification factors based on the water-level oscillations in the Well X10: (a) due to pore pressure fluctuations; (b) due to vertical displacement. A wellbore radius of rw = 54 mm was used in the calculation. Other parameters, such as storativity (S), porosity (n), wave velocity (v), etc., were all obtained using eqs (3) to (7). From the M7.8 Solomon Islands earthquake data, we obtained S = 19.2 × 10−6, n = 0.25, and hydraulic conductivity K = 0.375 cm s−1. 4 RESULTS AND DISCUSSIONS We can also acquire aquifer parameters based on the coseismic responses induced by nine other earthquakes (Fig. 1) using the methods described above. The results are shown in Table 1. Table 1. Estimated hydraulic parameters of Well X10 aquifer. Quantity  Unit  2016/12/25  2016/12/20  2016/12/17  2016/12/9  2016/12/8  2016/12/7  2016/11/25  2016/11/25  2016/11/22  2016/11/13  Magnitude  Mw  7.6  5.5  7.9  7.8  6.0  6.6  6.6  6.9  6.9  7.8  Distance  km  18,537  775  8451  9543  113  4348  1248  13,818  4520  12,798  Amplitude of WL  mm  15.0  5.0  4.0  27.5  30.0  3.5  4.5  24.5  18.5  15.0  Amplitude of VD  mm  0.582  0.057  0.109  1.046  0.430  0.094  0.168  0.604  0.488  0.582  Period  s  23  8  21  28  10  17  26  32  28  18  Wave velocity  m s−1  3461  3598  3732  3579  4302  3641  3821  3917  3725  3072  Wave Length  m  80 542  28 780  79 613  101 796  44 051  62 134  97 828  125 333  105 955  56 166  Shh  10−5  1923 983  6699  84 874  5495 827  324 297  98 186  454 044  158 855  1067 488  1944 934  Sww  10−5  1394  1  44  9813  186  26  262  270  1484  784  Shw  10−5  48,108  58  1853  223,039  6729  1478  10,003  6181  39,271  36,008  Storativitya  10−6  15.42  22.39  13.25  19.15  25.09  12.15  11.32  15.18  16.24  16.49  Storativityb  10−6  14.33  13.96  12.66  18.40  21.74  11.18  10.39  14.33  16.02  15.21  Storativityc  10−6  16.60  35.89  13.88  19.94  28.96  13.21  12.34  16.08  16.46  17.88  Storativityd  10−6  15.45  24.08  13.26  19.16  25.26  12.18  11.35  15.19  16.24  16.53  Specific Storage  10−6 m−1  1.93  3.01  1.66  2.40  3.16  1.52  1.42  1.90  2.03  2.07  Porositye  —  0.20  0.32  0.18  0.25  0.33  0.16  0.15  0.20  0.22  0.22  Hydraulic  cm s−1  0.05  0.10  0.09  0.38  0.09  0.09  0.07  0.10  0.12  0.08  conductivityf                        Transmissivity  m2 d−1  325  691  588  2592  622  608  449  677  802  518  Hydraulic  m2 s−1  243  332  513  1565  285  578  458  516  571  363  diffusivity  Quantity  Unit  2016/12/25  2016/12/20  2016/12/17  2016/12/9  2016/12/8  2016/12/7  2016/11/25  2016/11/25  2016/11/22  2016/11/13  Magnitude  Mw  7.6  5.5  7.9  7.8  6.0  6.6  6.6  6.9  6.9  7.8  Distance  km  18,537  775  8451  9543  113  4348  1248  13,818  4520  12,798  Amplitude of WL  mm  15.0  5.0  4.0  27.5  30.0  3.5  4.5  24.5  18.5  15.0  Amplitude of VD  mm  0.582  0.057  0.109  1.046  0.430  0.094  0.168  0.604  0.488  0.582  Period  s  23  8  21  28  10  17  26  32  28  18  Wave velocity  m s−1  3461  3598  3732  3579  4302  3641  3821  3917  3725  3072  Wave Length  m  80 542  28 780  79 613  101 796  44 051  62 134  97 828  125 333  105 955  56 166  Shh  10−5  1923 983  6699  84 874  5495 827  324 297  98 186  454 044  158 855  1067 488  1944 934  Sww  10−5  1394  1  44  9813  186  26  262  270  1484  784  Shw  10−5  48,108  58  1853  223,039  6729  1478  10,003  6181  39,271  36,008  Storativitya  10−6  15.42  22.39  13.25  19.15  25.09  12.15  11.32  15.18  16.24  16.49  Storativityb  10−6  14.33  13.96  12.66  18.40  21.74  11.18  10.39  14.33  16.02  15.21  Storativityc  10−6  16.60  35.89  13.88  19.94  28.96  13.21  12.34  16.08  16.46  17.88  Storativityd  10−6  15.45  24.08  13.26  19.16  25.26  12.18  11.35  15.19  16.24  16.53  Specific Storage  10−6 m−1  1.93  3.01  1.66  2.40  3.16  1.52  1.42  1.90  2.03  2.07  Porositye  —  0.20  0.32  0.18  0.25  0.33  0.16  0.15  0.20  0.22  0.22  Hydraulic  cm s−1  0.05  0.10  0.09  0.38  0.09  0.09  0.07  0.10  0.12  0.08  conductivityf                        Transmissivity  m2 d−1  325  691  588  2592  622  608  449  677  802  518  Hydraulic  m2 s−1  243  332  513  1565  285  578  458  516  571  363  diffusivity  aFrom Shh/Sww. bFrom Shh/Shw. cFrom Shw/Sww. daverage. eFrom Clark's method and eq. (7). fFrom Cooper's model. Notes: WL: water level; VD: vertical displacement. View Large Although the magnitudes and epicentres of the 10 earthquakes differ (the epicentral distances vary from 8451 to 18 537 km, which is almost half of the earth's perimeter), the estimated hydraulic parameters are similar. The porosity ranges from 0.15 to 0.33, the storativity ranges from 11.4 × 10−6 to 25.3 × 10−6 (corresponding to a specific storage Ss = S/d of 1.42–3.16 × 10−6 m−1), and the hydraulic conductivity ranges from 0.05 to 0.12 cm s−1 (corresponding to a transmissivity T = K·d of 300–800 m2 d−1). The hydraulic conductivity following the M7.8 earthquake on 2016 December 9 is 0.375 cm s−1, which is much larger than the other nine earthquakes. 4.1 Aquifer parameter changes induced by earthquakes Many researchers have noted that the aquifer permeability can be enhanced by seismic waves (Brodsky et al.2003; Elkhoury et al.2006; Manga et al.2012; Wang et al.2016), and this has also been proven by rock physics experiments (Elkhoury et al.2011; Candela et al.2014). Aquifer permeability enhancement always occurs with the coseismic water level abrupt/gradual response (Brodsky et al.2003; Sun et al.2015). The water-level changes in Well X10 were induced by the earthquakes from November—December 2016; however, nine of the ten showed oscillations but no subsequent step changes, indicating that the aquifer permeability did not change following the seismic waves (Shi et al.2015). A step-rise change occurred in the corrected water level following the M6.0 earthquake on 2016 December 8 (Fig. 8a), although the hydraulic conductivity did not increase immediately and showed an increasing water level after the earthquake. As a result, the estimated hydraulic conductivity during the M6.0 earthquake on 2016 December 8 was similar to the other earthquakes, and it showed substantial enhancement during the M7.8 earthquake on 2016 December 9, which was significantly higher than the other nine earthquakes. We speculated that the highest hydraulic conductivity estimated from the M7.8 earthquake data may be caused by the enhancement of hydraulic conductivity in the Well X10 aquifer after the M6.0 earthquake on 2016 December 8. Such enhancement, however, was caused by seismic wave shaking and new fractures induced by the stress/stain adjustment after the M6.0 earthquake (113 km from Well X10). The distribution of the coseismic strain field shows that Well X10 was in the compressive area of M6.0 earthquake and subjected to a coseismic strain of approximately 6 × 10−9 (Fig. 10), and the earthquake also has a higher seismic energy density than the other nine earthquakes (Fig. 3). As shown in Fig. 1a, Well X10 was drilled into a fault zone, which is especially sensitive to stress–strain variations. Figure 10. View largeDownload slide Coseismic strain field of the 2016 December 9 M7.8 Solomon Islands earthquake. Figure 10. View largeDownload slide Coseismic strain field of the 2016 December 9 M7.8 Solomon Islands earthquake. The coseismic strain was too small to increase the hydraulic conductivity, but the higher seismic energy density or the stress-strain adjustment after the M6.0 earthquake caused by the closer distance has the potential to create new fractures in the fault zone, which subsequently increases the porosity and hydraulic conductivity. As shown in Fig. 8(a), the step-rise recovery in the corrected water level (corresponding to the parameter changes) lasts for approximately 10 d and did not recover when the M7.8 earthquake occurred (12 hr after the M6.0 earthquake). 4.2 Model sensitivity to the estimated parameters The parameter estimations show that the storativity and hydraulic conductivity were less controlled by the earthquake magnitude and epicentral distance (Fig. 11, seismic energy density is closely related to the earthquake magnitude and epicentral distance) while the storativity was affected by the main seismic wave period. The storativity estimated from two earthquakes with smaller periods, the M6.0 earthquake on 2016 December 8 and M5.5 earthquake on 2016 December 20, was significantly greater than that for the other earthquakes, which may be related to the closeness between the significant seismic wave period (the corresponding period of the highest CSD at 10 s and 8 s of M6.0 and M5.5, respectively) and the resonant period of Well X10 $$(2\pi \sqrt {{H_e}/g} = 8.5\,{\rm{s}})$$. The amplification factor of the well-aquifer system will increase when the seismic wave period is closer to the well resonant period, which will increase the ratio of the power spectrum, Shh/Sww, Shh/Sww, and average storativity. However, the seismic wave period would not affect the hydraulic conductivity estimate because data from the period ranging from 5 to 50 s were all collected to estimate the hydraulic conductivity (Fig. 9), which was different from the estimation of storativity (only the significant period was selected). In addition, the highest estimated hydraulic conductivity in Fig. 11b on 2016 December 9 was caused by the enhancement of aquifer permeability after the M6.0 earthquake on 2016 December 8. Figure 11. View largeDownload slide The relation between seismic energy density and storativity (a), hydraulic conductivity (b), in which the colour circle shows the period of seismic wave. Figure 11. View largeDownload slide The relation between seismic energy density and storativity (a), hydraulic conductivity (b), in which the colour circle shows the period of seismic wave. Eqs (3)–(6) show that the storativity S was mainly affected by the wave length λ, aquifer thickness d and finite length L, where λ and d were determined by the observed data and L needs a manual setting. As shown in Fig. 12(a), different L values result in different curves of CSD, and the main difference is the smoothness of those curves. Although the smoothness was affected by the manual setting of L, the estimated results of porosity n and storativity S did not change significantly. Figure 12. View largeDownload slide Model sensitivity to the estimated parameters: (a) curve smoothness of CSD caused by different length (L) of data, and the corresponding porosity (n) and storativity (S); (b) estimated values of hydraulic conductivity (K) by inputting different wave velocity (v), porosity (n), wellbore radius (rw), storativity (S), or effective height of water column (He) in which each curve shows the changes of K caused by the increase of each input parameter while other parameter are unaltered. Figure 12. View largeDownload slide Model sensitivity to the estimated parameters: (a) curve smoothness of CSD caused by different length (L) of data, and the corresponding porosity (n) and storativity (S); (b) estimated values of hydraulic conductivity (K) by inputting different wave velocity (v), porosity (n), wellbore radius (rw), storativity (S), or effective height of water column (He) in which each curve shows the changes of K caused by the increase of each input parameter while other parameter are unaltered. From eqs (8) and (9), the amplification factors were determined by the wellbore radius rw, porosity n, wave velocity v, effective water volume height He (related to the thickness of the aquifer), storativity S and transmissivity T (or hydraulic conductivity K) if the other parameters were given, such as the bulk modulus of elasticity of water Ew, density of water ρ and gravitational acceleration g. As shown in Fig. 12(b), based on the data of the M7.8 earthquake on 2016 December 9 for Well X10, the estimated hydraulic conductivity using the model of Cooper et al. (1965) increases with an increase in rw, n and v and decreases with a decrease in S and He. Although the estimated hydraulic conductivity was affected by those input parameters, the magnitude of the estimated K did not change when the inputs were given as reasonable values. 4.3 Parameters estimated by different methods As shown in Fig. 8(a), the tidal response was clearer with the corrected water level than the raw water level, and a step-wise change induced by the M6.0 earthquake can also be identified. Many studies have employed the tidal response of the water level as a proxy for studying the hydraulic parameters of well-aquifer systems (Elkhoury et al.2006; Xue et al.2013; Shi & Wang 2016; Wang et al.2016), and the most commonly used method for this analysis is Hsieh's horizontal flow model (Hsieh et al.1987). When the aquifer is well confined, the permeability can be determined by the model using the negative phase shift between the water level and the tidal strain. However, Hsieh's model does not apply when the phase shift between the water level and the tidal strain is positive (Shi & Wang 2016). Roeloffs (1996) proposed that a positive phase shift results from the diffusion of pore pressure to the water table (the aquifer is not well confined), and the degree of confinement controls the frequency response of the pore pressure to tidal loading (Kinoshita et al.2015). Because hydraulic diffusivity is related to transmissivity and storativity, D = T/S, and the relationship between phase shift and hydraulic diffusivity can be written as follows:   $$p(z,\omega ) = B{K_u}\varepsilon (1 - {e^{ - (1 + i)z/\sqrt {2D/\omega } }})$$ (10)where p(z,ω) is the pore-pressure fluctuation at depth z, B is Skempton's coefficient, Ku is the bulk modulus of the saturated rock under undrained conditions, ɛ is the change in the volumetric strain and ω is the frequency of fluctuation. Roeloffs (1996) also mentioned that the response as a function of time (t) to a unit strain step (ɛ0) can be written as follows:   $$p(z,t) = - B{K_u}{\varepsilon _0}{\rm{erf}}(\sqrt {{z^2}/4Dt} ).$$ (11) By using the Venedikov's program of tidal analysis (Venedikov et al.2003), the phase shift of the water level in Well X10 was analysed. For the analysis, we set a data span of 28 d with a sliding step of one day. Then, we used the average value of the phase shifts of each month to determine the aquifer parameters. As shown in Fig. 13(a), the phase shift was obvious from October 2016 to January 2017. Based on Hsieh's horizontal flow model and Roeloffs’ vertical diffusion model, we can estimate the transmissivity and hydraulic diffusivity of aquifers (Shi & Wang 2016). To differentiate the parameters from the two models, we used subscript h to indicate the horizontal model and v for the vertical model. In addition, the water-level step-rise after M6.0 earthquake (marked by a dashed rectangle in Fig. 13a) was enlarged as shown in Fig. 13(b), and the hydraulic diffusivity was estimated according to eq. (11). Figure 13. View largeDownload slide Changes of water level and tidal phase shift in Xin 10 well during October 2016 to January 2017 (a). Water level and phase shift variations during October 2016 to January 2017 (b). Coseismic water level response following the 2016-12-08 M6.0 earthquake and fitting curve of the vertical model by using different hydraulic diffusivity values. Figure 13. View largeDownload slide Changes of water level and tidal phase shift in Xin 10 well during October 2016 to January 2017 (a). Water level and phase shift variations during October 2016 to January 2017 (b). Coseismic water level response following the 2016-12-08 M6.0 earthquake and fitting curve of the vertical model by using different hydraulic diffusivity values. The estimation results are listed in Table 2. Based on the coseismic response, the values of S and Th were estimated by the models mentioned in this study (Cooper et al.1965; Shih 2009); and based on the tidal response, the values of S and Th were estimated by Hsieh's horizontal flow model and Dh was calculated by Dh = Th/S. Then, based on Roeloffs’ vertical diffusion model, Dv was estimated by positive phase shifts and afterquake recovery and corresponding Tv was calculated by Tv = DvS. The results show that the value of storativity estimated from coseismic response was slightly smaller (less than an order of magnitude), although the value of transmissivity was much higher (approximately two orders of magnitude) than the value inferred from tidal response. The value of vertical diffusivity Dv estimated from tidal response was similar to the value from the afterquake recovery. Table 2. Aquifer parameters estimated from different data.     Transmissivity  Hydrulic diffusivity  Responses  Storativity (10−6)  Th (m2 d−1)  Tv (m2 d−1)  Dh (m2 d−1)  Dv (m2 d−1)  Coseismic response  10–20  320–800  —  2–5 × 107  —  Tidal response  70–120a  2–14a  0.002–0.004b  0.2–2 × 105a  30.83b  Stain step response  —  —  —  —  20–50      Transmissivity  Hydrulic diffusivity  Responses  Storativity (10−6)  Th (m2 d−1)  Tv (m2 d−1)  Dh (m2 d−1)  Dv (m2 d−1)  Coseismic response  10–20  320–800  —  2–5 × 107  —  Tidal response  70–120a  2–14a  0.002–0.004b  0.2–2 × 105a  30.83b  Stain step response  —  —  —  —  20–50  aFrom the negative phase shift in Fig. 13(a). bFrom the positive phase shift in Fig. 13(b). Notes: The em-dash sign (—) means it is not applicable. View Large The higher estimated value in transmissivity by the seismic wave method is not caused by an inappropriate model or error in estimation but rather to the different frequency response of the well-aquifer system. As with the barometric response (Fig. 8b), the response coefficient varies with the frequency of loading. For Well X10 aquifer system, it is likely that the higher the loading frequency, the greater the transmissivity (or hydraulic conductivity). Under periodic loading, especially in the short term, such as a seismic wave with a period less than one minute, the water in the aquifer might not flow far. As in Darcy's law (K = v/I), for a given hydraulic gradient I, periodic loading in the short term will speed up the flow of water (v), which will increase the effective hydraulic conductivity, K. Therefore, the transmissivity estimated from the coseismic response will be higher than the value from the tidal response. In addition to the response from coseismic and barometric loading, this mechanism is also reflected in the response to wind. From the water level (raw or corrected) shown in Fig. 8(a), fluctuations caused by short-term strong winds are also obvious. In addition, the parameter values estimated from the horizontal flow model were higher than the values from the vertical flow model (Table 2), which may indicate that the water flow in Well X10 aquifer is mainly horizontal. 4.4 Scale effect when estimating hydraulic parameters Accurate observations are important for obtaining reliable results. The Rayleigh wave timescale is approximately tens of seconds. Therefore, the sampling rate of the seconds scale is critical, followed by the resolution on the millimetre scale. As shown in Table 1, the water-level amplitudes induced by the seismic wave range from 3.5 to 30 mm, which requires the resolution of the pressure sensor be at least 1 mm. Moreover, the precision and accuracy of the sensor must be higher. In this study, the pressure sensors have a resolution of 1 mm and sampling rate of 1 s, which provide a guarantee for the complete record of water-level fluctuation. The sensor is located approximately 1 m below the surface of the water, with a precision of 0.02 per cent and accuracy of 0.05 per cent. The F.S. measurement error (repeatability of successive measurements under the same conditions) is approximately 0.2 mm, and the accuracy (the difference between the measured value and the true) is 0.5 mm. Because we are concerned with the relative changes in water level, the measurement error of 0.2 mm meets our research needs. To obtain accurate estimates of hydraulic conductivity, water must flow within the well-aquifer system. The water-level fluctuation in the well was caused by the periodic loading of seismic waves, and the motion of the water column created a periodic discharge from the aquifer into the well accompanied by drawdown, which is approximated by Cooper et al. (1965)   $$s = \frac{{\omega r_w^2{x_0}}}{{2T}}[Ker(\alpha )\cos (\omega t) - Kei(\alpha )\sin (\omega t)]$$ (9)where $$\alpha = r\sqrt {\omega S/T}$$ (dimensionless); x0 is the amplitude of the water-level fluctuation; r is the radial distance from the centre of well; and $$\omega r_w^2{x_0}/2T$$ represents the effect of the hydraulic gradient, which mainly affects the amplitude of drawdown. Ker(α)cos (ωt) − Kei(α)sin (ωt) represents the effective water storage of the aquifer; thus, it determined the range that produces significant drawdown. If we set S = 10−5, Tτ = 0.1 m2 and rw = 6 cm, the region of significant drawdown is approximated by a radius of 100 m. Hence, the parameters we estimate from coseismic water-level fluctuation can represent the hydraulic conductivity/transmissivity at a scale of hundreds of metres. In addition to the amplitude of water-level response to surface motion, the phase angle is as follows:   $$\eta = \arctan \left[\frac{{\omega r_w^2gKer({\alpha _w})}}{{2T{H_e}({\omega ^2} - \omega _w^2)}}\right]$$ (10)where $${\omega _w} \approx \sqrt {g/{H_e}}$$, which is the resonant frequency of a well. According to eq. (10), the phase ranges from −90° to 90° and decreases with an increase in frequency. The phase is negative when the frequency is less than the natural frequency and positive when the natural frequency is greater than the natural frequency; in addition, the phase is reversed when the frequency is equal to the natural frequency. As in the barometric efficiency/response method, the seismic wave and tide response methods that estimate hydraulic conductivity are also frequency dependent. Additionally, the hydraulic conductivity may vary over an order of magnitude on a small scale because it is dependent on the scale at which it is measured. Thus, discussing the scale effect on the method is necessary. Because the Cooper et al. (1965) and Hsieh et al. (1987) models were derived from the model of Theis (1935), their assumptions and solutions were similar. Here, if we assume that the amplitudes of periodic water-level fluctuation are the same, for example, 10 cm with a well radius of 6 cm, the amplitude of the discharge induced by the seismic wave will be on a dozens of cubic metres scale, whereas only dozens of litres for the tidal response would induce discharge. The former is more similar to the discharge produced in the field pumping tests. Because of the frequency scale, the model of Hsieh et al. is more suitable for determining aquifer transmissivities less than 1 m2 d−1; otherwise, the estimated value will be smaller than the actual value. The model of Cooper et al. is more suitable for aquifers with larger transmissivity. Hence, the transmissivity results, Th, inferred from the tidal response are smaller than the actual values as shown in Table 2. In addition, the large hydraulic diffusivity values of Dh may not represent the actual value of the whole aquifer system because of (1) the several hundred metres scale and periodic loading, and (2) the assumption that the release of groundwater from storage is instantaneous because of the drop in water head. 5 CONCLUSIONS This paper presents a method for estimating the hydraulic parameters of confined aquifers using water-level oscillations induced mainly by Rayleigh waves. The methodology is demonstrated using the coseismic responses of water levels in Well X10 to distant earthquakes and correlations with seismic waves that were analysed using spectral methods. The water-level oscillations in Well X10 are correlated with seismic waves in both their amplitude and frequency, and the aquifer hydraulic properties inferred from the seismic wave response of Well X10 yield a hydraulic conductivity of K = 1 × 10−3 m s−1 and a storativity of S = 15 × 10−6. The method proposed in this study is especially useful for wells that must be observed continuously and for which a clear tidal signal cannot be recorded. Acknowledgements The original water level and seismic wave data supporting the work presented in this paper are from the monitoring centre of the Earthquake Agency of the Xinjiang Uygur Autonomous Region. These data were extremely important to the study and the knowledge gained. We thank Todd C. Rasmussen and an anonymous reviewer for their constructive reviews and the help provided by the Editor. This research was supported by the National Natural Science Foundation of China (41502239 to Sun XL, U1602233 and 41602266 to Shi ZM), the National Key R&D Program of China (2017YFC1500502 to Sun XL) and the Research Grant from the Institute of Crustal Dynamics, China Earthquake Administration (ZDJ2017-27 to Sun XL). REFERENCES Allègre V., Brodsky E.E., Xue L., Nale S.M., Parker B.L., Cherry J.A., 2016. Using earth-tide induced water pressure changes to measure in situ permeability: a comparison with long-term pumping tests, Water Resour. Res. , 52( 4), 3113– 3126. https://doi.org/10.1002/2015WR017346 Google Scholar CrossRef Search ADS   Bendat J.S., Piersol A.G., 2000. Random Data, Analysis and Measurement Procedures , John Wiley and Sons. Brodsky E.E., Roeloffs E., Woodcock D., Gall I., Manga M., 2003. A mechanism for sustained groundwater pressure changes induced by distant earthquakes, J. geophys. Res. , 108( B8), doi:10.1029/2002JB002321. https://doi.org/10.1029/2002JB002321 Candela T., Brodsky E.E, Marone C., Elsworth D., 2014. Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing, Earth plant. Sci. Lett. , 392, 279– 291. https://doi.org/10.1016/j.epsl.2014.02.025 Google Scholar CrossRef Search ADS   Chapuis R.P., 1992. Similarity of internal stability criteria for granular soils, Can. Geotech. J. , 29( 4), 711– 713. https://doi.org/10.1139/t92-078 Google Scholar CrossRef Search ADS   Clark W.E., 1967. Computing the barometric efficiency of a well, J. Hydraul. Div. , 93( 4), 93– 98. Cooper H.H., Bredehoeft J.D., Papadopulos I.S., Bennett R.R., 1965. The response of well-aquifer systems to seismic waves, J. geophys. Res. , 70( 16), 3915– 3926. https://doi.org/10.1029/JZ070i016p03915 Google Scholar CrossRef Search ADS   Cooper H.H., Jacob C.E., 1946. A generalized graphical method for evaluating formation constants and summarizing well-field history, EOS, Trans. Am. geophys. Un. , 27( 4), 526– 534. https://doi.org/10.1029/TR027i004p00526 Google Scholar CrossRef Search ADS   Dawson K.J., Istok J.D., 1991. Aquifer Testing: Design and Analysis of Pumping and Slug Tests , CRC Press. Duhamel P., Vetterli M., 1990. Fast Fourier transforms: a tutorial review and a state of the art, Signal Process. , 19( 4), 259– 299. https://doi.org/10.1016/0165-1684(90)90158-U Google Scholar CrossRef Search ADS   Elkhoury J.E., Brodsky E.E., Agnew D.C., 2006. Seismic waves increase permeability, Nature , 441( 7097), 1135– 1138. https://doi.org/10.1038/nature04798 Google Scholar CrossRef Search ADS PubMed  Elkhoury J.E., Niemeijer A., Brodsky E.E., Marone C., 2011. Laboratory observations of permeability enhancement by fluid pressure oscillation of in situ fractured rock, J. geophys. Res. , 116( B2), doi:10.1029/2010JB007759. https://doi.org/10.1029/2010JB007759 Freeze R.A., Cherry J.A., 1979. Groundwater , Prentice-Hall. Hsieh P.A., Bredehoeft J.D., Farr J.M., 1987. Determination of aquifer transmissivity from earth tide analysis, Water Resour. Res. , 23( 10), 1824– 1832. https://doi.org/10.1029/WR023i010p01824 Google Scholar CrossRef Search ADS   Jacob C.E., 1940. On the flow of water in an elastic artesian aquifer, EOS, Trans. Am. geophys. Un. , 21( 2), 574– 586. https://doi.org/10.1029/TR021i002p00574 Google Scholar CrossRef Search ADS   Jiao J.J., Rushton K.R., 1995. Sensitivity of drawdown to parameters and its influence on parameter estimation for pumping tests in large-diameter wells, Ground Water , 33( 5), 794– 800. https://doi.org/10.1111/j.1745-6584.1995.tb00026.x Google Scholar CrossRef Search ADS   Kinoshita C., Kano Y., Ito H., 2015. Shallow crustal permeability enhancement in central Japan due to the 2011 Tohoku earthquake, Geophys. Res. Lett. , 42( 3), 773– 780. https://doi.org/10.1002/2014GL062792 Google Scholar CrossRef Search ADS   Ma L., Xu Y.S., Shen S.L., Sun W.J., 2014. Evaluation of the hydraulic conductivity of aquifers with piles, Hydrogeol. J. , 22( 2), 371– 382. https://doi.org/10.1007/s10040-013-1068-y Google Scholar CrossRef Search ADS   Manga M., Beresnev I., Brodsky E.E., Elkhoury J.E., Elsworth D., Ingebritsen S.E., Wang C.Y., 2012. Changes in permeability caused by transient stresses: field observations, experiments, and mechanisms, Rev. Geophys. , 50( 2), doi:10.1029/2011RG000382. https://doi.org/10.1029/2011RG000382 Rasmussen T.C., Crawford L.A. 1997. Identifying and removing barometric pressure effects in confined and unconfined aquifers, Groundwater , 35( 3), 502– 511. https://doi.org/10.1111/j.1745-6584.1997.tb00111.x Google Scholar CrossRef Search ADS   Roeloffs E., 1996. Poroelastic techniques in the study of earthquake-related hydrologic phenomena, Adv. Geophys. , 37, 135– 195. Google Scholar CrossRef Search ADS   Sethi R., 2011. A dual-well step drawdown method for the estimation of linear and non-linear flow parameters and wellbore skin factor in confined aquifer systems, J. Hydrol. , 400( 1), 187– 194. https://doi.org/10.1016/j.jhydrol.2011.01.043 Google Scholar CrossRef Search ADS   Shi Z., Wang G., 2016. Aquifers switched from confined to semiconfined by earthquakes, Geophys. Res. Lett. , 43( 21), 11 166–11 172. Google Scholar CrossRef Search ADS   Shi Z., Wang G., Manga M., Wang C., 2015. Mechanism of co-seismic water level change following four great earthquakes—insights from co-seismic responses throughout the Chinese mainland, Earth planet. Sci. Lett. , 430, 66– 74. https://doi.org/10.1016/j.epsl.2015.08.012 Google Scholar CrossRef Search ADS   Shih D.C.F., 2009. Storage in confined aquifer: spectral analysis of groundwater responses to seismic Rayleigh waves, J. Hydrol. , 374( 1), 83– 91. https://doi.org/10.1016/j.jhydrol.2009.06.002 Google Scholar CrossRef Search ADS   Shih D.C.F., Wu Y.M., Chang C.H., 2013. Significant coherence for groundwater and Rayleigh waves: evidence in spectral response of groundwater level in Taiwan using 2011 Tohoku earthquake, Japan, J. Hydrol. , 486, 57– 70. https://doi.org/10.1016/j.jhydrol.2013.01.013 Google Scholar CrossRef Search ADS   Stockwell R.G., Mansinha L., Lowe R.P., 1996. Localization of the complex spectrum: the S transform, IEEE Trans Signal Process , 44( 4), 998– 1001. Google Scholar CrossRef Search ADS   Sun X, Wang G, Yang X., 2015. Coseismic response of water level in Changping well, China, to the Mw 9.0 Tohoku earthquake, J. Hydrol. , 531, 1028– 1039. https://doi.org/10.1016/j.jhydrol.2015.11.005 Google Scholar CrossRef Search ADS   Theis C.V., 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, EOS, Trans. Am. geophys. Un. , 16( 2), 519– 524. https://doi.org/10.1029/TR016i002p00519 Google Scholar CrossRef Search ADS   Venedikov A.P., Arnoso J., Vieira R., 2003. VAV: a program for tidal data processing, Comput. Geosci. , 29 ( 4), 487– 502. https://doi.org/10.1016/S0098-3004(03)00019-0 Google Scholar CrossRef Search ADS   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( 8), 5923– 5937. https://doi.org/10.1002/2016WR018893 Google Scholar CrossRef Search ADS   Wang C.-Y., 2007. Liquefaction beyond the Near Field, Seismol. Res. Lett. , 78( 5), 512– 517. https://doi.org/10.1785/gssrl.78.5.512 Google Scholar CrossRef Search ADS   Wang C.-Y., Manga M., 2010. Hydrologic responses to earthquakes and a general metric, Geofluids. , 10( 1–2), 210– 216. Xue L., Li H.B., Brodsky E.E., Xu Z.Q., Kano Y., Wang H., Yang G., 2013. Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone, Science , 340( 6140), 1555– 1559. https://doi.org/10.1126/science.1237237 Google Scholar CrossRef Search ADS PubMed  © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial