## Abstract

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 M_{2} and O_{1} 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 M_{2} and the O_{1} 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.

## 1. Introduction

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 M_{2} wave and O_{1} 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. Data and Methods

### 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 (*M _{w} 7.9*) in May 2008 and the Tohoku earthquake (

*M*) 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.

_{w}9.1### 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

A tidal analysis was performed with the M_{2} wave (0.5175 day) and the O_{1} 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 M_{2} wave and the O_{1} 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\u2032$. If there was no leakage in the aquitard ($K\u2032=0$), the mixed flow model would be equivalent to Hsieh’s horizontal flow model [10].

_{2}and O

_{1}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 M

_{2}was divided into that of O

_{1}to omit $BKu$; therefore, three inputs were used to obtain three parameters ($S$, $T$, and $K\u2032$):

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

## 3. Results

### 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\xd710\u22126$ to $9.0\xd710\u22126$ m^{2}/s, indicating the horizontal permeability of the aquifer increased. The storage coefficient of the aquifer $S$ increased from $8.0\xd710\u22126$ to $3.5\xd710\u22125$ without significant changes. The vertical hydraulic conductivity of the aquitard $K\u2032$ remained at $1.0\xd710\u221210$, indicating that no vertical leakage occurred in the aquitard. In contrast, after the Tohoku earthquake, the transmissivity increased from $1.8\xd710\u22126$ to $1.0\xd710\u22125$ m^{2}/s, the storage decreased from $1.5\xd710\u22125$ to $1.0\xd710\u22125$, and the vertical hydraulic conductivity increased from $1.0\xd710\u221210$ to $2.0\xd710\u22127$, 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.

Hydraulic parameters | Wenchuan | Tohoku | ||||
---|---|---|---|---|---|---|

Preearthquake | Postearthquake | Standard error | Preearthquake | Postearthquake | Standard error | |

Horizontal transmissivity ($T$) | 2.20$E$-06 | 9.00$E$-06 | 2.52$E$-06 | 1.80$E$-06 | 1.00$E$-05 | 7.63$E$-07 |

Storage coefficient ($S$) | 8.00$E$-06 | 3.50$E$-05 | 9.37$E$-06 | 1.50$E$-05 | 1.00$E$-05 | 3.78$E$-06 |

Vertical hydraulic conductivity ($K\u2032$) | 1.00$E$-10 | 1.00$E$-10 | 0 | 1.00$E$-10 | 2.00$E$-07 | 5.76$E$-08 |

Hydraulic parameters | Wenchuan | Tohoku | ||||
---|---|---|---|---|---|---|

Preearthquake | Postearthquake | Standard error | Preearthquake | Postearthquake | Standard error | |

Horizontal transmissivity ($T$) | 2.20$E$-06 | 9.00$E$-06 | 2.52$E$-06 | 1.80$E$-06 | 1.00$E$-05 | 7.63$E$-07 |

Storage coefficient ($S$) | 8.00$E$-06 | 3.50$E$-05 | 9.37$E$-06 | 1.50$E$-05 | 1.00$E$-05 | 3.78$E$-06 |

Vertical hydraulic conductivity ($K\u2032$) | 1.00$E$-10 | 1.00$E$-10 | 0 | 1.00$E$-10 | 2.00$E$-07 | 5.76$E$-08 |

$\u2217$Standard errors are calculated with the data of 4 months before the earthquakes for $T$ and $S$, while with the data of 4 months after the Tohoku earthquake was used for $K\u2032$.

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\xd710\u22126$ to $1.8\xd710\u22126$ m^{2}/s, and the vertical hydraulic conductivity $K\u2032$ remained at less than $1\xd710\u221210$ 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. Discussion

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

According to the above equation, we know that the phase shift is simultaneously affected by three parameters—$S$, $T$, and $K\u2032$. Figure 6(a) is the solution space of the M_{2} wave phase shift of the mixed flow model. The same phase shift corresponds to multiple groups of $S$, $T$, and $K\u2032$ 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\xd710\u22126$ m^{2}/s, $S=1\xd710\u22125$, and $K\u2032=1\xd710\u22127$ to $1\xd710\u22126$ m/s, whereas the M_{2} wave phase shift was still negative. In addition, when the $K\u2032$ was very small ($K\u2032$ ≈0, approximately equal to less than $1\xd710\u221210$ 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\u2032$ value was less than $1\xd710\u22128$ m/s, with a specific value of the storage coefficient $S$, the phase shift would not be affected even if $K\u2032$ continued to decrease. This phenomenon indicates that there may be a threshold for $K\u2032$.

Figure 6(b) is the curve of $K\u2032$ and the phase shifts of the M_{2} and O_{1} waves. When$\u2009K\u2032$ was small enough (≤$1\xd710\u22128$ m/s), the phase shifts almost remained unchanged. Thus, $K\u2032=1\xd710\u22128$ m*/s* may be the threshold value indicating no obvious vertical flow occurred.

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

### 4.2. Relationship of the Phase Shifts of the M_{2} and O_{1} Waves for the Mixed Flow Model

For the mixed flow model with the fixed reasonable values of $S$, $T$, and $K\u2032$, the tidal phase shifts will vary with the tidal frequency. Figure 7 demonstrates the relationship of the phase shifts of the M_{2} and O_{1} waves with the same parameters. The phase shift of the M_{2} wave should theoretically always be less than that of the O_{1} 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 M_{2} wave adjusted to be equal to or larger than that of the O_{1} wave during the time spans studied.

The probable explanations for this are as follows:

- (1)
The O

_{1}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\u2032$ are all within a reasonable range. The calculated results of $K\u2032$ 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 M_{2} 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\xd710\u22125$ 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\u2032$ 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 M_{2} 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 O_{1} and M_{2} 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\u2032$) 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\u2032$. 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 $Ss$$Ss=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} m^{2}, the transmissivity $T$ can be calculated form the permeability using the relationship the transmissivity $T=b\u2217\rho gk/\mu $, where $\rho $ and $\mu $ are the water density (1000 kg/m^{3}) and viscosity (0.001 Pa$\u2217$s) [10, 26]. The reference range of $T$ is about 10^{-5}-10^{-3} m^{2}/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\u2032$) 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 *M _{w}* 7.9 earthquake and 2011 Tohoku

*M*9.1 earthquake in Zuojiazhuang area and found that the Tohoku vertical vibration energy density was greater than 0.1 J/m

_{w}^{3}; meanwhile, the Wenchuan seismic energy density was less than 0.1 J/m

^{3}. 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 O

_{1}phase shift is greater than the M_{2}phase shift, this model can accurately calculate the parameters. When two phase shifts are similar, or the M_{2}phase shift is greater than the O_{1}phase shift, there may be errors in the calculation. However, it can still roughly estimate the order of magnitude of hydraulic parameters

## 5. Conclusions

As both the M_{2} and O_{1} 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\u2032$ value increased from $1.0\xd710\u221210$ to $2.0\xd710\u22127$, 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.

## Data Availability

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

## Conflicts of Interest

The authors declare that they have no conflicts of interest.

## Acknowledgments

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).

*Reviews of Geophysics*

*Water Resources Research*

*Water Resources Research*

*Chinese Science Bulletin*

*Earth and Planetary Science Letters*

*Tidal analysis of borehole pressure-A tutorial*

*Journal of Geophysical Research: Solid Earth*

*Geology*

*Water Resources Research*

*Journal of Geophysical Research: Solid Earth*

*Journal of Volcanology and Geothermal Research*

*Journal of Geophysical Research: Solid Earth*

*Geophysical Research Letters*

*Water Resources Research*

*Geophysical Research Letters*

*Science*

*Water Resources Research*

*Distribution of Active Faults in China (1:4000000)*

*Seismology Analysis and Prediction System Based on GIS (Mapsis Software)*

*Journal of Geophysical Research: Solid Earth*

*Journal of Geophysical Research: Solid Earth*

*Water Resources Research*

*Theory of linear poroelasticity with applications to geomechanics and hydrogeology*

*Characterizing small-scale permeability of the Arbuckle Group, Oklahoma*

*Rockmass Mechanics (in Chinese)*

*Journal of Seismology*