Permeability is an important parameter of shallow crust that affects underground water migration. It is important to calculate the permeability changes caused by earthquakes for underground engineering safety and pollutant dispersion. The response analysis of groundwater in wells, to Earth tides, has been widely used to calculate permeability. However, the calculation results and the judgment of leakage are controversial. In this study, the water level response of both the M2 and O1 tidal waves in the Zuojiazhuang well was analyzed. An optimization fitting, using the mixed flow model, was performed to calculate the hydraulic parameters of aquifers and aquitards. After the Wenchuan earthquake, the horizontal transmissivity of aquifers increased temporarily and then slowly recovered. However, after the Tohoku earthquake, the horizontal transmissivity increased permanently by one order of magnitude. The vertical hydraulic conductivity increased permanently by three orders of magnitude, meaning that vertical leakage occurred. As both the M2 and the O1 waves are used as constraints for the calculation with the mixed flow model, the results of this study are more precise than those of previous ones. This article also gives some suggestions for utilizing the mixed flow model.

Permeability is a key point in underground engineering. Research shows that earthquakes can cause changes in the permeability of the crust [1–3]. According to Hsieh’s horizontal flow model, the horizontal transmissivity of aquifers can be estimated from water level records of Earth’s solid tides [1, 4–6]. For aquitards, vertical hydraulic conductivity can be calculated by the vertical flow model [3, 7–9]. A general model for the tidal response of groundwater in aquifers, presenting both horizontal flow and vertical leakage, was established on the basis of these two idealized models [10].

Earthquakes can cause water level changes (e.g., [11, 12]). The main reason why earthquakes can change water levels in aquifers is that the dynamic stress of seismic waves can clear the sediment inside the cracks and lead to an increase in permeability (e.g., [1, 13]). There are also multiple examples of permeability reduction (e.g., [14, 15]). Additionally, the permeability in some regions can slowly recover after earthquakes (e.g., [16, 17]), whereas some regions may be permanently changed (e.g., [6]). For the Zuojiazhuang (ZJZ) well, Zhang et al. [6] calculated the hydraulic parameters of the aquifer by the radial flow model and the vertical flow model. Zhang et al. [18] used the barometric response function method to calculate the permeability in both the aquifer and the aquitard, and the mixed flow model was used to assist in the identification of the vertical hydraulic conductivity.

In this study, the mixed flow model was used independently to calculate the hydraulic parameters in the aquifer and aquitard by the tidal responses of water levels in the ZJZ well. First, a tidal analysis was performed to extract the M2 wave and O1 wave from the water level data and Earth tide data. Next, a tidal response analysis was completed to analyze the changes in the amplitude ratio and phase shift. Finally, optimization was conducted by fitting the mixed flow model and the obtained hydraulic parameters. The results show that the vertical hydraulic conductivity increased by three orders of magnitude after the 2011 Mw 9.1 Tohoku earthquake. It is necessary to consider the impact of vertical leakage in future underground engineering, so that the safety of workers could be guaranteed and the leakage of the underground resources and environment pollutions could be monitored.

2.1. Observation

The ZJZ well is located in Beijing, China (116.45° E, 39.95° N) (Figure 1(a)) and is near the Shunyi-Liangxiang fault [18, 19] with a depth of 2,605 m. The radius of the casing is 73 mm, and the well hole is 80 mm. Figure 1(b) is the lithologic well logging of the ZJZ well. This study researched the influence the Wenchuan earthquake (Mw 7.9) in May 2008 and the Tohoku earthquake (Mw 9.1) in March 2011 had on the water level of the ZJZ well because both caused phase shifts and steep changes in the water levels. Both were far-field earthquakes, with an epicentral distance of 1,354.68 km for the Wenchuan earthquake and 1,532.22 km for the Tohoku earthquake.

2.2. Data

For the ZJZ well, groundwater occurs in the rock fractures. Water level observations for the ZJZ well began in late 2007. The water level was recorded by the LN-3A water level gage produced by the Seismological Research Institute of the China Earthquake Administration with an accuracy of ≤0.2% full scale, a resolution of 1 mm, and a sampling rate of 1 per minute.

As shown in Figure 2, the period selected was six years, from January 2008 to January 2015. Water level data are not valid after 2015 owing to silt cleaning in the well. The Earth tide data utilized were from Mapsis, a software package that has been tested against other well-known software packages such as Baytap-G [20].

2.3. Tidal Analysis

Earth tides have many frequency components. The response is the sum of all the frequencies:
(1)yth=∑k=wavesak∗eiωkt,
where ak is the complex coefficients. The response of the well is
(2)yresp=∑k=waveshk∗ak∗eiωkt.
A tidal analysis was performed to extract the complex coefficients hk in each frequency, including amplitudes and phases. According to Doan’s research, assuming that hk was the same for frequencies ωk close to each other, the equation could be simplified as
(3)yresp=∑g=groupshg∗∑k=wavesak∗eiωkt,
where hg is a complex parameter and a key point in the tidal analysis [7]. The results of the tidal analysis were amplitude ratios and phase shifts between water levels and tides, caused by the elastic properties of rocks and the time it took for water to enter and discharge from the pores of those rocks.

A tidal analysis was performed with the M2 wave (0.5175 day) and the O1 wave (1.0758 days) waves. There are two reasons these waves were analyzed. First, their amplitudes were large enough to be easily extracted. Second, they were created only by gravity and not affected by barometric loads [4, 7, 21, 22]. The Fourier transform was performed to obtain the frequency series of the water level and Earth’s tides. Once the corresponding components in the frequency series were obtained, the inverse Fourier transform was performed to return to the time series. Note that the frequency of the interference wave is near the M2 wave and the O1 wave in the frequency series. The length of the time series should be considered in component extraction (see Appendix A of [6]). A moving time window of 30 days and a running step of 3 days were used for the continuous calculation.

2.4. Model

The mixed flow model (Figure 3) was used to simulate the system of the aquifer and aquitard. This model has three independent parameters: the storage coefficient S⁠, transmissivity T⁠, and vertical hydraulic conductivity K′⁠. If there was no leakage in the aquitard (⁠K′=0⁠), the mixed flow model would be equivalent to Hsieh’s horizontal flow model [10].

Using the mixed flow model, the amplitude ratio A and phase shift η can be defined as
(4)A=abshw,0/BKuε0ρg,(5)η=arghw,0/BKuε0ρg,
or
(6)A=absiωSiωS+K′/b′ξ,(7)η=argiωSiωS+K′/b′ξ,
here
(8)ξ=1+rcrw2iωrw2TβK0βrwK1βrw,(9)β=K′Tb′+iωST1/2,
where hw,0 is the water level tidal wave, ε0 is the theoretical solid-earth tidal strain wave, A is the relative amplitude ratio of the water level to solid-earth tidal strain, η is the phase shift, arg and abs are the angle and magnitude of complex numbers, respectively, ω is the circular frequency of the tidal wave, BKu is the product of Skempton’s coefficient and the bulk modulus, Kn is the second kind of modified Bessel function (BesselK function) with order n⁠, K′ is the vertical hydraulic conductivity of the aquitard, T is the horizontal transmissivity of the aquifer, S is the storage coefficient of the aquifer, b′ is the thickness of aquitard, rc is the case radius, and rw is the well radius.
The absolute amplitude ratio A′=hw,0/ε0 was directly obtained from the water level and solid-earth tides, so the equation set can be rewritten as
(10)A′=absiωSiωS+K′/b′ξ∗BKuρg ,η=argiωSiωS+K′/b′ξ,
where S⁠, T⁠, K′⁠, and BKu are the unknown parameters. The original equation set had four unknowns and two equations. Both the M2 and O1 tidal responses (phase shifts and amplitude ratios) were used to be the four inputs of those equations, since the specific value of BKu cannot be confirmed. Instead, the amplitude ratio of M2 was divided into that of O1 to omit BKu⁠; therefore, three inputs were used to obtain three parameters (⁠S⁠, T⁠, and K′⁠):
(11)A′M2/AO1′=absiωM2SiωM2S+K′/b′ξ/absiωO1SiωO1S+K′/b′ξ,ηM2=argiωM2SiωM2S+K′/b′ξ,ηO1=argiωO1SiωO1S+K′/b′ξ.

The least squares optimal fitting method was used to resolve those equations, and the most optimal solution was obtained.

3.1. Water Level Response to Earth Tides

The tidal analysis was performed with the time window of 30 days and the running step of 3 days. As shown in Figure 4, the phase shifts and amplitude ratios increased after the earthquakes. The phase shifts were negative before and after the Wenchuan earthquake but became positive after the Tohoku earthquake. Leakage cannot be determined only by positive or negative phase shifts [10]. A conclusion on leakage can be reached after the calculation of the hydraulic parameters.

3.2. Hydraulic Parameters Obtained from the Mixed Flow Model

As shown in Figure 5, for the Wenchuan earthquake, the horizontal transmissivity T increased from 2.2×10−6 to 9.0×10−6 m2/s, indicating the horizontal permeability of the aquifer increased. The storage coefficient of the aquifer S increased from 8.0×10−6 to 3.5×10−5 without significant changes. The vertical hydraulic conductivity of the aquitard K′ remained at 1.0×10−10⁠, indicating that no vertical leakage occurred in the aquitard. In contrast, after the Tohoku earthquake, the transmissivity increased from 1.8×10−6 to 1.0×10−5 m2/s, the storage decreased from 1.5×10−5 to 1.0×10−5⁠, and the vertical hydraulic conductivity increased from 1.0×10−10 to 2.0×10−7⁠, suggesting that significant vertical leakage occurred in the aquifer. All of these changes remained permanent over time. Table 1 shows the comparison of hydraulic parameters before and after the earthquakes.

Normally, the permeability of an aquifer increases after seismic waves pass through, likely owing to the unclogging of the fractures induced by the flush of the fast flow fluids [1, 11]. If the energy of the seismic wave did not reach the threshold of a permanent change in the stratigraphic structures, the permeability will recover slowly over time. After the Wenchuan earthquake and during the period between these two earthquakes, the horizontal transmissivity T slowly recovered from 9.0×10−6 to 1.8×10−6 m2/s, and the vertical hydraulic conductivity K′ remained at less than 1×10−10 m/s. After the Tohoku earthquake, the hydraulic parameters did not return to preearthquake levels, which meant the seismic waves had changed the physical properties of the stratigraphic structures permanently. This is similar to the conclusion reached by Zhang et al. [6].

4.1. Optimal Solutions and the Threshold of the Occurrence of Vertical Flow

Nonlinear functions of equation (11) have multiple solutions.
(12)ηM2=argiωM2SiωM2S+K′/b′ξ.

According to the above equation, we know that the phase shift is simultaneously affected by three parameters—S⁠, T⁠, and K′⁠. Figure 6(a) is the solution space of the M2 wave phase shift of the mixed flow model. The same phase shift corresponds to multiple groups of S⁠, T⁠, and K′ values. Therefore, it is necessary to select and judge the best optimal solution with the geological data (value range of the parameters of the lithology of the aquifer and aquitard). From the solution space, the traditional viewpoint of judging the leakage of the aquifer by the positive or negative phase shifts is defective. For example, obvious vertical flow happened with T=1×10−6 m2/s, S=1×10−5⁠, and K′=1×10−7 to 1×10−6 m/s, whereas the M2 wave phase shift was still negative. In addition, when the K′ was very small (⁠K′ ≈0, approximately equal to less than 1×10−10 m/s), the mixed flow model could be approximated to the horizontal flow model [10].

In the solution space of Figure 6(a), it is easy to find that when the K′ value was less than 1×10−8 m/s, with a specific value of the storage coefficient S⁠, the phase shift would not be affected even if K′ continued to decrease. This phenomenon indicates that there may be a threshold for K′⁠.

Figure 6(b) is the curve of K′ and the phase shifts of the M2 and O1 waves. When K′ was small enough (≤1×10−8 m/s), the phase shifts almost remained unchanged. Thus, K′=1×10−8 m/s may be the threshold value indicating no obvious vertical flow occurred.

A comparison was needed to verify whether 1×10−8 m/s was the threshold. Figure 6(c) is the hydraulic parameter solutions of the mixed flow model when the K′ values were fixed at 0, 1×10−10⁠, 1×10−9⁠, and 1×10−8 m/s. For the horizontal transmissivity T⁠, when K′ values were less than 1×10−8 m/s, the solutions were nearly the same. But for the storage coefficient S⁠, only when K′ is equal to or less than 1×10−10 m/s were the solutions consistent. Considering both solutions of S and T⁠, the K′ value of 1×10−10 m/s could be used as the threshold of the occurrence of vertical flow in the aquitard. Therefore, when K′ is equal to or less than 1×10−10 m/s, no vertical flow happens.

4.2. Relationship of the Phase Shifts of the M2 and O1 Waves for the Mixed Flow Model

According to equation (7),
(13)ηx=argiωxSiωxS+K′/b′ξ.

For the mixed flow model with the fixed reasonable values of S⁠, T⁠, and K′⁠, the tidal phase shifts will vary with the tidal frequency. Figure 7 demonstrates the relationship of the phase shifts of the M2 and O1 waves with the same parameters. The phase shift of the M2 wave should theoretically always be less than that of the O1 wave.

As shown in Figure 4(a), these relationships were broken after the Wenchuan earthquake and the Tohoku earthquake, and the phase shifts of the M2 wave adjusted to be equal to or larger than that of the O1 wave during the time spans studied.

The probable explanations for this are as follows:

  • (1)

    The O1 wave is unstable and has relatively large errors, which could be observed from the tidal analysis results

  • (2)

    The mixed flow model is an idealized model and has many simplifications and assumptions [10]. These simplifications and assumptions are as follows: the aquifer and the aquitard are assumed to be homogeneous and isotropic media; the aquifer is laterally extensive, whereas the flow in the aquitard is vertical; the aquifer should be significantly thinner than the aquitard; and the aquitard is assumed to be incompressible and has zero storage. Therefore, many complex in situ geological factors have not been considered in the mixed flow model

  • (3)

    From the most related analytical solution of the mixed flow model of the ZJZ well (Figure 8), it can be seen that the phase shifts and the corresponding results of S⁠, T⁠, and K′ are all within a reasonable range. The calculated results of K′ are in the same order of the expected magnitude. Therefore, the results could roughly estimate the order of magnitude of the vertical hydraulic conductivity of the aquitard and judge whether there is any vertical leakage. The results may not be absolutely accurate, but they are still acceptable and instructive

4.3. Comparison with Previous Results

The ZJZ well has been researched previously. This study’s results have been compared to the results of Zhang et al. [6] and Zhang et al. [18] and other available data for reference.

Zhang et al. [6] used the tidal response method to calculate the hydraulic parameters of the aquifer of the ZJZ well from 2008 to 2015 with the vertical flow model and the horizontal flow model, respectively. Before the Tohoku earthquake, the phase shift of the M2 wave was maintained as negative. They used the horizontal flow model to do the calculation. After the Tohoku earthquake, the phase shift turned positive or near zero. These results showed that the seismic waves of the Tohoku earthquake had permanently changed the physical properties of the strata medium of the ZJZ well. This research used traditional ideal models and could be improved upon by using the mixed flow model.

Zhang et al. [18] actually utilized a mixture of models and methods to calculate the hydraulic parameters from 2009 to 2012. Firstly, they used the barometric response model to calculate the static-confined barometric efficiency (⁠BE⁠), the horizontal transmissivity of aquifer (⁠Taqu⁠), the vertical hydraulic diffusivity of this confining layer (⁠Dcon⁠), and the vertical pneumatic diffusivity of this zone (⁠Dunsat⁠); however, the storativity of the aquifer (⁠Saqu⁠) is determined by the tidal analysis with the pure vertical flow model and explained obscured in the paper (and the storativity of the confining layer (⁠Scon⁠) is roughly used 1×10−5 in the calculation, which has little influence on the results). Next, they put these obtained parameters (⁠Taqu and Saqu⁠) into the mixed flow model of Wang et al. [10] and achieved the vertical conductivity of the aquitard (⁠K′ value). This value must be in accordance with the results of the previous barometric responses since they used the same input parameters (⁠Taqu and Saqu⁠) and as indicated by Wang and Manga [23] the analytical solutions of the mixed flow model [10] and that of the barometric response model for the mixed flows [24] are the same.

Additionally, both studies only used the M2 wave to do the tidal analysis. Their results will not be as constrained as this study’s results where the mixed flow model was utilized with both the O1 and M2 waves.

The mixed flow model [10] is the latest superlative model that can be used in the research on tidal response of the layer, which has the same analytical solution as the barometric response model of Rojstaczer [24] for mixed flows [23]. Previous studies could only use the horizontal flow model [4] or the vertical flow model [7] to calculate the hydraulic parameters of the layer. The mixed flow model has many simplifications and assumptions, including 2-layer model and considering the entire overlayer as an aquitard, and the vertical hydraulic conductivity (⁠K′⁠) reflects the average permeability of the entire overlayer. In order to verify the rationality of our results, we refer to some existing rock physical test data, when calculating the optimal solution of S⁠, T⁠, and K′⁠. As provided in Figure 1(b), the borehole logs are available. Both the dolomite and limestone are carbonate rocks. Referring to Wang [25], the typical value for the specific storage SsSs=S/b is about 10-7-10-5 m-1. Our calculation result is that storage coefficient S is about 10-5-10-4 and the aquifer thickness is 521 m that means Ss is about 10-7-10-6 m-1. Meanwhile, the permeability of dolomites k is about 10-14-10-12 m2, the transmissivity T can be calculated form the permeability using the relationship the transmissivity T=b∗ρgk/μ⁠, where ρ and μ are the water density (1000 kg/m3) and viscosity (0.001 Pa∗s) [10, 26]. The reference range of T is about 10-5-10-3 m2/s. Both S (and Ss⁠) and T (and k⁠) which we calculated are reasonable values. In addition, according to the results of Liu and Tang [27], the hydraulic conductivity K of dense andesite and limestone with no cracks is about between 10-10 and 10-8 m/s, and our result of vertical hydraulic conductivity (⁠K′⁠) of the overlayer is also within this reasonable range before the occurrence of the 2011 Tohoku earthquake.

Our research shows that, after the Tohoku earthquake, the hydraulic parameters did not return to preearthquake levels, which meant the seismic waves had changed the physical properties of the stratigraphic structures permanently. Zhang et al. [6] calculated the seismic wave energy of 2008 Wenchuan Mw 7.9 earthquake and 2011 Tohoku Mw 9.1 earthquake in Zuojiazhuang area and found that the Tohoku vertical vibration energy density was greater than 0.1 J/m3; meanwhile, the Wenchuan seismic energy density was less than 0.1 J/m3. Only for the energy of Tohoku earthquake exceeded the energy density threshold for permanent changes in the layer [6]. Therefore, it is highly possible to cause permanent changes in the physical properties of the shallow surface layer media in the Zuojiazhuang area, as well as the permeability. However, because of the 2-layer simplified assumptions of the mixed flow model which perhaps is the latest superlative mixed flow model, so far, the overlayer (aquitard) is seemed as a whole layer, and we could not know which layer within the aquitard was damaged by seismic waves detailly. In addition, the coseismic tensile stress and reopening of preexisting fractures may also be one of the reasons for the hydraulic parameters change [28]. In the future study, we may research on how earthquakes promote the development of fractures and the permeability with more wells and testing.

4.4. Suggestions for Using the Mixed Flow Model

  • (1)

    Most strata are calculated with unknown BKu⁠; thus, three input parameters with three output result parameters are proper and reasonable

  • (2)

    Multiple solutions exist in the calculation results; thus, the initial testing value range for the least squares fitting is important, as it is needed to set the lithology of the aquifers and aquitards as the constraints

  • (3)

    In the case where the O1 phase shift is greater than the M2 phase shift, this model can accurately calculate the parameters. When two phase shifts are similar, or the M2 phase shift is greater than the O1 phase shift, there may be errors in the calculation. However, it can still roughly estimate the order of magnitude of hydraulic parameters

As both the M2 and O1 waves are used as constraints for the calculation with the mixed flow model, these results are more precise than those from previous works. After the Wenchuan earthquake in May 2008, no vertical leakage occurred in the aquitard of the ZJZ well. The horizontal permeability of the aquifer increased after the earthquake and then slowly recovered and decreased to preearthquake levels. However, after the Tohoku earthquake in 2011, the K′ value increased from 1.0×10−10 to 2.0×10−7⁠, indicating that vertical leakage occurred in the aquitard and the horizontal permeability of the aquifer increased permanently. Using the mixed flow model is a reasonable way to calculate the hydraulic parameters. Compared with traditional models, it can calculate the vertical hydraulic parameters of the aquitard more accurately and can precisely judge whether vertical leakage occurred. However, because of the simplifications and ideal assumptions of the mixed flow model [10], it may be different from the in situ situation, and there may be calculation errors.

Groundwater level data used in this study were from the Groundwater Monitoring Network of the China Earthquake Networks Center (CENC): doi:10.12197/2021GA006.

The authors declare that they have no conflicts of interest.

Groundwater level data used in this study were from the Groundwater Monitoring Network of the China Earthquake Networks Center (CENC): doi:10.12197/2021GA006. We thank Professor Chi-yuen Wang, Rui Yan, Shengle Li, Wang Zhang, Xiqiang Zhou, and Yuchuan Ma for helpful discussions. Research was supported by the National Natural Science Foundation of China (41874161), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant Number XDA14010303), and the Youth Innovation Promotion Association Foundation of the Chinese Academy of Sciences (2019069).