Possible Continuous Vertical Water Leakage of Deep Aquifer: Records from a Deep Well in Tianjin Province, North China

It is well known that the storage of the valuable helium and other natural gases needs thick well-con ﬁ ned aquitards. Meanwhile, the assumption that deep aquifers are con ﬁ ned is most evident from the wide practices of targeting deep aquifers as the storage of toxic wastes and CO 2 . With the negative phase shifts of the M 2 tidal wave, previously, deep Gaocun well ( ~ 3500m) in the North China Platform is assumed to be con ﬁ ned fairly well. However, water level of the Gaocun well has been continuously decreasing for ~ 40 years without coseismic variations and without explainable mechanism. In this study, innovatively, we ﬁ nd that even buried in ~ 3500 meters deep, and covered by thick compact mudstones, the aquifer of well Gaocun is calculated to be continuously leaking for ~ 40 years probably induced by the vertical leakage incurred by aquitard fractures, which is under the assumptions of the none-leaking bottom layer and none direct pumping or exploiting in the observation aquifer layer. Meanwhile, for the ﬁ rst time, we did a detailed systematic comparison between the leaky models of the tidal response and the barometric response of water level, which indicate a consistency leaky result, and the minor di ﬀ erence mainly induced by the frequency di ﬀ erences. Merits and demerits of both models are also analyzed. Last but not least, underground ﬂ uid leakages might frequently occur, which are becoming increasingly worldwide urgent since that might also be related with surface sourced contaminations, toxic waste burials, and burial and exploration of the natural gas reservoirs.


Introduction
Deep aquifers are always assumed to be confined, which is most evident from the worldwide practices of targeting deep aquifers as the storage of toxic wastes (e.g., [1][2][3]) and for the proposed storage of sequestrated CO 2 (e.g., [4]). Moreover, Benson and Cole [4] explicitly stated that aquifers deeper than 800 m and confined by thick and extensive aquitard should be suitable for the sequestration of CO 2 , whereas, nowadays, groundwater leakage is an incremental vital problem for many critical aquifers worldwide (e.g., [1]), and safety of the underground waste repositories might be influenced by variations of vertical permeability in confined layers. Plus, the agricultural waste, such as the pesticides and fertilizers, also might migrate from the surface soil down into the groundwater through the vertical leakage fractures. Hence, sustained observation of the confinement of aquifers should be needed to detect any leakage to allow appropriate prevention of the contamination of groundwater resource on the one hand and the outflow of the stored wastes into the environment on the other.
Traditionally, pure horizontal flow model [5] or pure vertical flow model [6,7] is commonly used to estimate hydraulic or poroelastic parameters. However, many aquifers may neither be perfectly confined nor fully unconfined but behave somewhere between the two endmembers [8].
For the tidal response of borehole water level to a semiconfined aquifer, Wang et al. [9] obtained a general analytical solution for the tidal response of water level with both horizontal flow and vertical flow (mixed-flow model). That is deduced from the early analytical solution of Hantush and Jacob [10], which is for the leaky and multilayered aquifers. Recently, the derived solutions were applied to study pressure variation in response to water injection/extraction in boreholes (e.g., [11]) or used to calculate hydraulic parameters both for the aquifer and for the aquitard (e.g., [9,[12][13][14]). Meanwhile, for the barometric response of borehole water level in a semiconfined aquifer, several analytical models have been deduced, and the analytical solution of Rojstaczer [15] was widely applied to estimate hydraulic properties (e.g., [13,[15][16][17]). In addition, Odling et al. [18] compared the numerical and analytical calculated results and also indicated that properties of a layer in an area extending some hundreds of meters from the monitored borehole could be estimated. Finally, as indicated by Wang and Manga [19], for the leaky aquifer, both the tidal response leaky model of Wang et al. [9] and the barometric response model of Rojstaczer [15] are with the same analytical solution, and the recent numerical testing of Zhu and Wang [20] indicates some limitations of this ideal analytical solution.
Although both the tidal response and the barometric response of borehole water level have been studied previously, the results obtained from the two models separately have never been compared quantitatively and deeply before. Though Zhang et al. [13] compared the results from both models, they only used the M 2 tidal wave for analysis, and they preset the obtained parameters (transmissivity and storativity of the aquifer) from the barometric response model as the input premise parameters for the tidal response leaky model, which must obtain the similar result of the vertical leakage (K ′ m/s) since as indicated by Wang and Manga [19] both of the two models are with the same analytical solution. In other words, they combined the barometric and the tidal response models together for the calculation, and the process is not independent. While for the study of Sun and Xiang [17] or Sun et al. [21], the obtained result from the barometric pressure analysis based on the model of Rojstaczer [15] has not calculated the coherence coefficient between the water level and the barometric pressure, and thus, the wide span (10 -2 cpd~102 cpd) of fitting should be unreasonable for all those wells since the coherence coefficients are very low (towards 0) when extended to the frequency region of >~10 cpd for most wells.
In this study, in order to test the confinement of the deep aquifer and to compare the tidal response and the barometric response leaky models independently, well Gaocun is set as the study objective. We used those two response (barometric response and tidal response) models totally separately to do calculations. Results obtained from the leaky model of the tidal response [9] solely with both the O 1 (diurnal) and M 2 (semidiurnal) tidal waves and that obtained from the barometric response model [15] are compared, which indicate similar and consistent leaky results, and the small differences for the value of the vertical leakage might be induced by different frequencies. Therefore, we justified that even at a depth of~3500 m, although the aquifer is covered with thick compact mudstone layers, significant leakage is occurring sustainably for~40 years. For the first time, deep, quantitative, and systematic comparisons between the leaky model of the tidal response and the barometric response leaky model have been made in this paper.

Methodology
2.1. Well Selection and Observation. The water-level data documented for the deep Gaocun well (3402.8 m) ( Figure 1 and Table 1) in the North China Platform (Figure 2(a)) indicates a continuously decreasing pattern without any coseismic variations for a long time (~40 years) even meeting with huge earthquakes (2008 M w 7.9 Wenchuan and 2011 M w 9.1 Tohoku earthquakes) (Figure 1). There are no reasonable explanations for this pattern, since the well was previously assumed to be confined fairly well because the M 2 tidal responses of the water level always had negative phase shifts ( Figure 3). Besides, the mudstone (clay) layer exists in the aquitard over the aquifer, which perhaps confirm the aquifer to be confined well. Hence, identifying the true mechanism was meaningful and significant. Well Gaocun is with the distance of greater than 100 km away from the ocean (Figure 2(a)); thus, oceanic tides could be avoided, which exert an impaction tens of kilometers away from the coast [22]. In addition, this well is in a relatively stable plain (e.g., the North China Platform), and thus, the outside disturbances and interferences are limited. The lithologic well log is provided in Figure 2 The basic information of this well is given in Table 1. From the year 1979, the SW40-1 (Chong-Qing) analog observation water-level instrument was used, which has a sampling rate of once per hour and a resolution of 1 mm.
Since the year 2001, the LN-3A digital water-level instrument was applied to document the water level, with an observation accuracy of ≤0.2% full scale, a resolution of 1 mm, and a sampling rate of once per minute. While the barometric pressure is documented with the RTP-1 digital instrument, the sampling rate is 1 sample per minute. For well Gaocun, the barometric observation starts from the year 2015. The software Mapsis was used [24] to calculate and analyze the theoretical Earth tides, and it has been tested against BAYTAP-G (e.g., [14]). In order to be uniform, we used the hourly data of water level from January 1st, 1981, to May 9th, 2020, and resampled into the minute data to do calculations in this study. For the barometric analysis, we used data from the year 2017 to the year 2019.

Tidal Response Analyses.
In this study, we used both O 1 and M 2 tide responses in the water-level data as constraints, and these have relatively large amplitudes and are free of barometric effects (e.g., [5,6,25]). For the analysis of the O 1 and M 2 phase shifts of the tidal strain and water level, a moving time window of 15 days (the appendix) and running step of 30 days are applied. Different window sizes have been tested, and the results are very similar. Since the time serious is~40 years long, the dots could be too many and gather together; thus, the relatively long running step (30 days) was applied. The observed M 2 and O 1 phases had opposite signs for the Gaocun well ( Figure 3). A negative phase response occurred with the M 2 wave, whereas a positive phase response occurred with the O 1 wave ( Figure 3). Strangely, the coseismic phase shifts (also water level ( Figure 1)) of this well remained unchanged ( Figure 3) even with huge earthquakes (the 2008 M w 7.9 Wenchuan earthquake or 2011 M w 9.1 Tohoku earthquake). Figure 4 shows the coherence estimation between the water level and the barometric pressure. It shows that the coherence was relatively high at low frequencies (< 0.8 cpd, vertical red line) while it deteriorates at high frequencies (> 0.8 cpd) because of the tidal and noises disturbances. Therefore, we concentrate on 0.01 cpd~0.8 cpd frequency range in the below analysis. The barometric response of groundwater is obtained by cross-spectrum analysis of the water-level and the barometricpressure time series. The time series of water-level was preprocessed to remove the effects of precipitation, tides, and recharge before the analysis [16].

Barometric Response Analyses.
The detailed analysis of the transfer functions is the same as the study of Zhang et al. [26]. The barometric pressure and water level were split in spans of 2N samples, so as to reduce noise error. Similar to the dealing of Zhang et al. [13], to be in accordance with the results of Rojstaczer [15] and Sun and Xiang [17], we reduced the barometric phase shifts by 180°, which is just a custom of Rojstaczer [15].

Parameters Obtained from the Tidal and
Barometric Response Leaky Models 3.1. Parameters Obtained from the Leaky Model with the Tidal Responses of the Water Level. Equations of the leaky model of Wang et al. [9] are shown as follows: Here,

Geofluids
A ′ is the absolute amplitude ratio between the water level and solid-earth tides, η is the phase shift, arg and abs are the angle and magnitude of complex numbers, ω is the circular frequency of the tidal wave, 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, r c is the case radius, and r w is the well radius. The subscripts M 2 and O 1 indicate the issues corresponding to the M 2 /O 1 tidal frequencies.
Hence, 3 input parameters (A M2 ′ /A O1 ′ , η M2 , and η O1 ) turn out 3 output parameters (S, T, and K ′ ), which will be more stable [14]. Since the least squares fitting relies on the input initial testing values of T, S, and K ′ , the aquifer lithology (each lithol-ogy has value spans for the hydraulic parameters [27]) and the aquitard lithology (reflect the vertical leakage extent K ′) are set as limitation conditions during the fitting processes. The detailed calculation processes are also introduced in Zhang et al. [14] and Yang et al. [12]. Briefly, during the non-linear, least squares inversion (Fsolve, https://www.mathworks.com/ help/optim/ug/fsolve.html), the initial values are set as the constrain conditions for the results, thus the obtained S and T are constrained by the parameter ranges of the aquifer rocks, while K' is the vertical leakage of the fracture, thus need to be >=0 as the constrain, and as the leakance term, K' should be relatively small according to the model setting of Wang et al. [9].
The analytical solution of the leaky model for well Gaocun is given in Figure 5. Finally, we obtained three parameters from 5 Geofluids the leaky model: the storage coefficient S and transmissivity T of the aquifer and the vertical conductivity K ′ of the aquitard ( Figure 6). From Figure 6, we can see the vertical leakage occurred even since 1981, without any occurrence of earthquakes, and the leakage sustained for~40 years, with the relatively stable vertical hydraulic conductivity (in accordance with the vertical permeability). Meanwhile, the horizontal transmissivity T (in accordance with the horizontal permeability) also shows a relatively stable value range but with a very small increase after the year 2001, which perhaps induced by unclogging of fracture clays ( Figure 6). There is no trade-off problem of the parameters in the inversion, and the inversion is independent and without preassumptions. Detailed information and the process of solutions are explained in the previous paper [14].
Zhang et al. [14] have done a deeply comparison between the calculated value based on the tidal response leaky model and the laboratory value for 6 wells on the North China Platform, including well GC, which show the reasonability of those obtained parameters.

Parameters
Obtained from the Barometric Response Leaky Model. As described in Sun and Xiang [17], the model of the barometric response of water level [15,16] also could be used to calculate numerous parameters of the aquifer and aquitard, and thus, the vertical leakage parameters also could be testified and to do comparisons with the above parameters obtained from the leaky model.
Similar as the data processed by Sun and Xiang [17], we also used the time series of water level and barometric pressure of January 1st, 2017, to May 30th, 2019, to do analysis. As we analyzed, only model C [17] could be strictly applied to well Gaocun (Figure 7), indicating a semiconfined aquifer, which is different from the fitting of Sun and Xiang [17] since lower frequency domain fitting (0.01 cpd~0.8 cpd) has been used in this study and is much more reasonable according to the coherence coefficient ( Figure 4).
Before the calculation of the barometric response, the storage coefficient of aquifer S aqu and the storage coefficient of aquitard S con need to be set in advance. The same as the dealing of Sun and Xiang [17], we also assigned S con = 10 −4 [13,15] and S aqu = nγβL aqu /BE [17,28] Table 2; we set n = 0:2 (Table 2) in this study [17,29], which is confirmed by the aquifer lithology (limestone), while γ is the weight density of water (9:8 × 10 3 N/m 3 ), β is the compressibility of water (4:8 × 10 −10 m 2 /N), and L aqu is the thickness of the aquifer, respectively. During the least squares fitting, the initial values are set as the constrain conditions for the results, and thus, the obtained T is constrained by the parameter range of the aquifer rock, while K ′ is the vertical leakage of the fracture, thus needs to be ≥ 0 as the constraint (e.g., Sun and Xiang [17]). Therefore, the result of the vertical leakage of the aquitard K ′ and the horizontal transmissivity of the aquifer T obtained from the barometric response leaky model are close to the result obtained from the leaky model based on the tidal responses of water level ( Table 2; Section 3.1).

Opposite Signs for the Phase Shifts of the M 2 and O 1 Tides Indicating Existing of Fractures
As studied by Hanson and Owen [30], the fracture orientation technique based on the tidal response of water level is an effective and convenient methodology which has no limitation of depth, and the fracture solutions of both natural and manually stimulated fractures can be obtained with static observation wells.
Traditionally, negative phase shifts indicate confined behavior (e.g., [5]) while positive phase shifts indicate unconfined behavior (e.g., [31]), but this is an oversimplification. The negative M 2 wave phase shifts led several researchers and engineers in the seismic prediction department of the China Earthquake Networks Center (CENC) to assume that the aquifer of the Gaocun well is confined fairly well. However, our findings indicated a very different situation. The opposite signs for the M 2 and O 1 responses (Figures 3 and 5) indicate substantial leakage (K ′ : 10 -7~1 0 -6 (m/s)) ( Figure 5) in the analytical solution of the tidal response leaky model. Meanwhile, as studied by Hanson and Owen [30], for a given vertical fracture, the O 1 and M 2 tides also display phase shifts with opposite signs. Similarly, the opposite signs for the phase shifts of the M 2 and O 1 tides can be explained by the fracture model of Bower [32], where the M 2 and O 1 responses are complementary with respect to the fracture strike and always have opposite signs under undrained conditions [32]. These support the robustness of the present conclusion that the inclined fracture exists in the aquitard of well Gaocun.
However, since the fracture model of Hanson and Owen [30] indicates a dry fracture without fluid or permeability, thus it is definitely different with the leaky model and could not be used to calculate the fault/fracture solutions in this paper.

Comparison between the Tidal/Barometric
Response Leaky Models  Figure 7: Fitting the (a) actual gain (barometric efficiency)-ω and (b) phase-ω curves calculated from the observed data of water level and barometric pressure with model C [17]. * In order to be in accordance with the result of Rojstaczer [15], the same as done by Zhang et al. [13], we reduced 180 degrees in the barometric phase shift. 7 Geofluids model with the vertical flow in the confining layer and the horizontal flow in the aquifer (Figure 8). The analytical solution of both models [9,15] has inherent specific similar simplifications: the assumption of the isotropic and homogeneous aquifer and aquitard separately; the aquitard storage was negligible which indicates a time constant for the aquitard to attain hydraulic equilibrium short compared to the period of the tidal or the barometric forcing; none basement leakage, and thus, the bedrock should be seemed as more stable and impermeable compared to the other layers; periodic discharge of the well will not cause impact on the water table, and the well is a line source (for the leaky model [9], the pore pressure in the well is the same with different radii but with the same depth, thus line source).
However, in fact, some wells could not strictly meet with the rigorous assumptions and only could do rough estimations.

Same Analytical Solutions of Those 2 Models.
For the best fitted model C, which equals to the aquifer with a semiconfined aquitard and neglecting the effect of the unsaturated zone [15,17], the differential equation of the barometric response is obtained by Rojstaczer [15].
As summarized by Wang and Manga [19], for the semiconfined aquifer, the differential equation of the barometric response model (model C of Rojstaczer [15]) is the same as that of the tidal response leaky model [9]. In fact, the analytical solution of both the two models is derived from the solution of Hantush and Jacob [10], which is for aquifer response to pumping in the leakance condition, and the difference is the well discharges at a periodic rate and without a constant rate.

Different Sensing Frequency Domains for Those 2
Models. However, some conditions are very different for those 2 response models, which perhaps induced the minor differences of the vertical hydraulic conductivity K ′ of the confining layer and the horizontal transmissivity T of the aquifer (Tables 3 and 4).
For groundwater response to tidal strain, we mainly used the diurnal (O 1 ) (~0.93 cpd) and semidiurnal (M 2 ) (~1.93 cpd) tidal waves to do the calculation, while the barometric response model is mainly based on the fitting with the relatively low-frequency domain dots, since the coherence coefficient is only relatively high in 0.03 cpd~0.8 cpd (Figure 4), and thus, the sensing frequency is very different, which perhaps induced the different results of calculations.

Merits and Disadvantages of Those 2 Models
. From the leaky model with tidal response, a sequence of S, T, and K ′ varied with time could be obtained, and with the moving calculation time windows, the minimum continuous data of only 15 days is needed. While based on the barometric response model, only one mean value of those parameters (T and K ′) could be obtained from the best fitting during a long-time span (at least~60-100-day continuous data to obtain one stable fitting value).
The value of the storativity of the aquifer S aqu and aquitard S con should be preset before the fitting and calculation with the barometric response model; thus, this calculation is not completely independent, and very probably the value of S aqu could not be determined precisely since it might vary with time, although which would not cause obvious impaction on the results of the other parameters (e.g., [16]). Accordingly, the inverted specific storage based on the tidal response leaky model of Wang et al. [9] (S aqu~1 0 -6 ( Table 3)) is much smaller than the laboratory obtained value S aqu = 10 −4 to 10 -3 (S s = 10 −6 for limestone and L aqu = 717:3 m for well Gaocun) [27]. The inconformity perhaps is induced by simplifications of the analytical model for the inversion [14], although this method is definitely independent or perhaps induced by the lack of the laboratory testing of more S values for different rocks, since it is not so important as the values of permeability. Therefore, the model of Wang et al. [9] is more flexible (only need observation of water level and need not to set the observation of the barometric pressure), independent (need not to do value preset for the parameters such as S aqu and S con ), and time continuity (could obtain timevaried sequences of parameters and could reflect the variations of parameters in relatively short time span (minimum time span~15 days for one value)). However, this tidal response leaky model calls for the relatively small aquifer thickness (thickness of the aquifer should be smaller than that of the aquitard), and the vertical flow in the aquitard should be smaller than the horizontal flow in the aquifer thus could be treated as the leakance term [10].
Thus, both the 2 models have their merits and demerits and with the same analytical solution.  Figure 8: Sketch showing the groundwater flow in response to varied barometric pressure, modified from Wang and Manga [19]. K ′ is the vertical hydraulic conductivity of the aquitard, and K is the hydraulic conductivity of the aquifer. 8 Geofluids

Discussion
In layered aquifers, both phase lag and attenuation data might be used to estimate hydraulic parameters of the aquifer, and inconsistent results may be due to propagation bias [33]. However, both the tidal response and the barometric response leaky models are with only one layer of aquifer (homogeneous and isotropic), and thus, the layered aquifer model is definitely different, hence could not be used to do comparisons.
6.1. Impaction Factors for the Groundwater Decrease. Numerous of incidents could cause the water-level decrease of well Gaocun, such as groundwater exploitation-related pumping in the shallower layer and precipitation decrease. However, as shown in Figure 9, as recorded by Tianjin Water Resource Bulletin, the precipitation amount of Tianjin City does not show a decrease pattern and thus could be excluded as the reason of the water-level decrease. Furthermore, previous studies indicate that groundwater depletion in North China Plain (NCP) is due to massive use of groundwater-based irrigation and the industry usage, and groundwater pumping caused the fall in the groundwater level and land subsidence (e.g., [34,35]). The NCP includes one shallow unconfined aquifer (40-60 m) and three deep confined aquifers of different depths (120-170 m, 250-350 m, and 400-600 m) [35], and the groundwater exploitation in the North China is always above the depth of 400 m [36]. Hence, the aquifer of Gaocun well is much deeper (3400 m) than the pumping layers (<400 m); thus, pumping of the groundwater in shallow layers might induce the deep groundwater of well Gaocun (3400 m) to migrate up to the shallower layer through the vertical fractures in the aquitard induced by the high pore pressure of the deeper layer and so as to reach a new balanced state.
In addition, the study of Zhao et al. [37] indicates that the subsidence rates change severely in time and space ( Figure 10) even in the same area of Tianjin. Therefore, the efforts to do the subsidence or groundwater exploitation statistics for the whole aquifer layer of well Gaocun or in the nearby megacities would be unrealistic and meaningless.

Reasonability of the Leaky Aquifer Model Encountering
Pumping in the Shallower Layer. Firstly, the pumping/ exploitation of the groundwater in the shallow layer primarily might induce the replenishment of the horizontal flow in the shallower layer and then might also induce the vertical migration of the deep groundwater, which could be very minor compared to the horizontal replenishment. Secondly, for the calculation of the tidal response leaky model, we distilled the O 1 (diurnal) and the M 2 (semidiurnal) tidal waves, through which the information of groundwater exploitation has been deleted already, except for those incidents with the same frequencies as the M 2 and O 1 waves, which rarely occurs. Therefore, we could still use the leaky model to deal with the leaky aquifer calculations under the circumstances of shallower layer pumping incidents in NCP.

Direct Pumping or Exploitation in the Observation
Aquifer Layer or Downward Leakage Induced by Fractures of the Bedrock Is Not Considered Here. The GC well (aquifer depth: 2685.5 m~3402.8 m) is close to the BD well (aquifer depth: 210 m~427 m) in the horizontal projection (37.6 km apart). Figure 11 shows the original water level of well GC and well BD. Both the 2 wells are with vertical leakage as calculated by Zhang et al. [14], with the K ′ ≈ 10 −5 (m/s) for well BD and K ′ ≈ 10 −7 (m/s) for well GC. Obviously, the water level decreased in both the 2 wells but decreased more slowly in well BD than in well GC, which probably induced by the horizontal fluid supplement for the aquifer of well BD in the shallow layer by rivers or lakes. Since there is hardly supplement of the fluid in the~3500 m deep aquifer, therefore, with the inclined-vertical fractures, the water level of well GC decreased more quickly than well BD.
Although it is well known that the Cambrian and Ordovician aquifer is also used for geothermal exploiting, as we investigated, there is no direct hydrothermal exploitation in the layer of the observation aquifer of well GC; as the government required, all those adjacent private hydrothermal exploitations also have been forbidden because of the ground subduction (Tianjin Geothermal exploration, Development and Design Institute); besides, most persons might think that with the  9 Geofluids gravity, the water should leak downward; however, as in the research region of hydrology, there is the basic rock (most likely granite and those tough rock) as the bottom layer of the aquifer for each well, which could be seemed to be without leakage. Thus, with the high temperature and high pore pressure in the deep aquifer, leaking only assumed to occur upward to the above confining layer. Therefore, direct pumping or exploitation in the observation aquifer layer or downward leakage induced by fractures of the bedrocks is not considered here in both the 2 models (tidal response leaky model and barometric response leaky model), which will make things much complex and different. 11 Geofluids (Ideal Similar Assumptions of Those 2 Models). Those ideal simplifications might induce some inversions to be unstable, such as kick points in Figure 6 and biased fitting in Figure 7(b).
Numerical simulation developed by Zhu and Wang [20] indicated that leakage of the basement and aquitard storage has significant impacts on the tidal response of borehole water level and the neglection of those factors in the analytical solution [9] could have undervalued the aquifer leakage.

Conclusion
For the first time, we did a detailed systematic comparison between the leaky models of the tidal response and the barometric response of water level, which indicates a consistency leaky result. Besides, we analyzed the merits and demerits of both the 2 models. With the confirmation of the existing of fractures and especially with the obtained K ′ calculated from both the tidal and the barometric response models, we verify that well Gaocun has a sustained obvious vertical leakage in the aquitard for~40 years, under the assumptions of the none-leaking bottom layer and none direct pumping or exploiting in the observation aquifer layer. The M 2 and O 1 waves need to be distilled with the same moving time window to facilitate comparison over the same time span. Spectral leakage between the M 2 (~0.0013419 circle/min) and S 2 (~0.0013889 circle/min) tides and between the O 1 (~0.0006455 circle/min) and K 1 (~0.0006963 circle/ min) (S 1 ) tides poses challenges for data analysis. Generally, the minimum time window separating the frequencies of the M 2 and S 2 tides in the frequency spectrum (i.e., two separate peaks) is set to be 29.5036 days (e.g., [26,38]). Meanwhile, 15 days are also proved to be suitable and reasonable (Appendix A in [26]). However, the frequencies of the O 1 and K 1 (S 1 ) tides cannot be separated over the same time window of 29.5036 days [26]. Similar to Appendix A of Zhang et al. [26], we found that both the O 1 and M 2 waves could be distilled precisely with a moving time window of 15 days (or 30 days). Time series for water level in wells always have transitions, missing data, or interruptions ( Figure 1). Thus, continuous waterlevel data with the same trend do not cover a sufficiently long time. A moving time window of 15 days is a relatively short time span that could reflect variations in the permeability more precisely.
Separation of the O 1 and K 1 tidal waves is as follows.
Tidal response analysis requires a frequency division: where df is the frequency difference, f s is the sampling frequency, N is the sampling points, and t is the sampling time [39]. The sampling frequency f s is 1 circle/min for the water level and Earth's solid tides in this paper.
To separate the O 1 and K 1 waves in the frequency spectrum (i.e., distinguish two wave peaks) (Figure 12(a)) and avoid spectral leakage (i.e., adjacent waves should be recorded integrally), df should be half the frequency difference df = ð 0:0006963 − 0:0006455Þ/2 = 2:54 × 10 −5 circle/min. In addition, the corresponding time window should be greater than or equal to T = 1/df = 27:34 days. The greatest common divisor for the frequencies of the O 1 and K 1 waves also needs to be calculated. With a calculated frequency of ≤2:54 × 10 −5 circle/ min, the accurate period should be 27.95 (27. and 15 days (Figure 14(b)) based on the data of Gaocun well before and after the 2008 Wenchuan earthquake. The results are similar. Thus, similar to the distillation of the M 2 wave, the 15-day time window can also be used to distill the O 1 wave.

Data Availability
The groundwater level data used in the analysis were from the Groundwater Monitoring Network of the China Earthquake Networks Center (CENC), which can be accessed from WDC for Geophysics, Beijing (http://wdc.geophys.ac .cn, doi:10.12197/2021GA001).

Conflicts of Interest
The authors declare that they have no conflicts of interest.