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

Authors: Yan Zhang; Xiaolong Sun; Tianming Huang; Shengwen Qi; Li-Yun Fu; Qiu-Ye Yang; Junhua Hu; Bowen Zheng; Wang Zhang
Journal: Geofluids
 (2022)

---

## Abstract

It is well known that the storage of the valuable helium and other natural gases needs thick well-confined aquitards. Meanwhile, the assumption that deep aquifers are confined is most evident from the wide practices of targeting deep aquifers as the storage of toxic wastes and CO2. With the negative phase shifts of the M2 tidal wave, previously, deep Gaocun well (~3500 m) in the North China Platform is assumed to be confined 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 find 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 first 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 difference mainly induced by the frequency differences. Merits and demerits of both models are also analyzed. Last but not least, underground fluid 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.


---

## Body

## 1. 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–3]) and for the proposed storage of sequestrated CO2 (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 CO2, 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–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–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 M2 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 O1 (diurnal) and M2 (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.


## 2. Methodology

### 2.1. Well Selection and Observation

The water-level data documented for the deep Gaocun well (3402.8 m) (Figure1 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 Mw 7.9 Wenchuan and 2011 Mw 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 M2 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(b).Figure 1
Original water level of well Gaocun (from January 1st, 1981, to May 9th, 2020) with an obvious downward trend. The gaps indicate invalid original data records. The vertical dashed lines represent the start time of the 2008Mw 7.9 Wenchuan earthquake and the 2011 Mw 9.1 Tohoku earthquake, which had not caused any influence on the water level. (GC is a nonartesian well; thus, the y-coordinate axis indicates the distance from the water-level surface to the ground.)Table 1
Basic information of well Gaocun.
Station name (abbreviation)
Epicentral distance (km)
Well diameter (mm)
Well depth (m)
Open interval (m)
Aquifer thickness (m)
Aquitard thickness (m)
Major lithology
Wenchuan
Tohoku
Outer
Inner
Gaocun (GC)
1478.00
1894.02
240.00
180.00
3 402.81
2685.50~3402.80
717.30
2685.50
LimestoneFigure 2
(a) Location of the North China Platform, the position of the groundwater well Gaocun (green triangle, GC: Gaocun), and the epicenters of the 2008Mw 7.9 Wenchuan earthquake and the 2011 Mw 9.1 Tohoku earthquake (from Global CMT). The “beach balls” show the focal mechanism for the 2 huge earthquakes. Red lines indicate the surface faults [23]. (b) Well logs of Gaocun. The red lines and numbers indicate the depth of the observation aquifer (the screened section). For analyzing the impaction of major rocks in the observation aquifer, we neglected those rocks with small amount (references: Seismic Monitoring Records of China; Seismological Press, Beijing, 2004 (in Chinese)).
(a)
(b)Figure 3
The phase responses over time (from January 1st, 1981, to May 9th, 2020) for well Gaocun at the M2 wave (red dots) and O1 wave (black dots) frequencies. The horizontal dashed line represents the 0° phase shifts, from which we can see the obviously separated M2 and O1 phase shifts. The vertical dashed lines represent the start time of the 2008 Wenchuan and the 2011 Tohoku huge earthquakes, which had not caused obvious influence on the phase shifts or the amplitude ratios. ∗We calculated the errors based on the data from Jan. 1st, 1981, to Dec. 31st, 1990 (which is relatively stable and continuous), and obtained the error of the phase shifts of the M2 wave with 4.08°, while O1 wave with 4.32°, and obtained the error of the amplitude ratios of the M2 wave with 65400.07 (m/strain), while O1 wave with 75700.15 (m/strain).The basic information of this well is given in Table1. 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.


### 2.2. Tidal Response Analyses

In this study, we used both O1 and M2 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 O1 and M2 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 M2 and O1 phases had opposite signs for the Gaocun well (Figure 3). A negative phase response occurred with the M2 wave, whereas a positive phase response occurred with the O1 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 Mw 7.9 Wenchuan earthquake or 2011 Mw 9.1 Tohoku earthquake).


### 2.3. Barometric Response Analyses

Figure4 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.Figure 4
Coherence coefficient of the barometric pressure and the processed water level. Data series of January 1st, 2017, to December 31th, 2017, were used to do coherence analysis.The barometric response of groundwater is obtained by cross-spectrum analysis of the water-level and the barometric-pressure time series. The time series of water-level was preprocessed to remove the effects of precipitation, tides, and recharge before the analysis [16].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].


## 2.1. Well Selection and Observation

The water-level data documented for the deep Gaocun well (3402.8 m) (Figure1 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 Mw 7.9 Wenchuan and 2011 Mw 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 M2 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(b).Figure 1
Original water level of well Gaocun (from January 1st, 1981, to May 9th, 2020) with an obvious downward trend. The gaps indicate invalid original data records. The vertical dashed lines represent the start time of the 2008Mw 7.9 Wenchuan earthquake and the 2011 Mw 9.1 Tohoku earthquake, which had not caused any influence on the water level. (GC is a nonartesian well; thus, the y-coordinate axis indicates the distance from the water-level surface to the ground.)Table 1
Basic information of well Gaocun.
Station name (abbreviation)
Epicentral distance (km)
Well diameter (mm)
Well depth (m)
Open interval (m)
Aquifer thickness (m)
Aquitard thickness (m)
Major lithology
Wenchuan
Tohoku
Outer
Inner
Gaocun (GC)
1478.00
1894.02
240.00
180.00
3 402.81
2685.50~3402.80
717.30
2685.50
LimestoneFigure 2
(a) Location of the North China Platform, the position of the groundwater well Gaocun (green triangle, GC: Gaocun), and the epicenters of the 2008Mw 7.9 Wenchuan earthquake and the 2011 Mw 9.1 Tohoku earthquake (from Global CMT). The “beach balls” show the focal mechanism for the 2 huge earthquakes. Red lines indicate the surface faults [23]. (b) Well logs of Gaocun. The red lines and numbers indicate the depth of the observation aquifer (the screened section). For analyzing the impaction of major rocks in the observation aquifer, we neglected those rocks with small amount (references: Seismic Monitoring Records of China; Seismological Press, Beijing, 2004 (in Chinese)).
(a)
(b)Figure 3
The phase responses over time (from January 1st, 1981, to May 9th, 2020) for well Gaocun at the M2 wave (red dots) and O1 wave (black dots) frequencies. The horizontal dashed line represents the 0° phase shifts, from which we can see the obviously separated M2 and O1 phase shifts. The vertical dashed lines represent the start time of the 2008 Wenchuan and the 2011 Tohoku huge earthquakes, which had not caused obvious influence on the phase shifts or the amplitude ratios. ∗We calculated the errors based on the data from Jan. 1st, 1981, to Dec. 31st, 1990 (which is relatively stable and continuous), and obtained the error of the phase shifts of the M2 wave with 4.08°, while O1 wave with 4.32°, and obtained the error of the amplitude ratios of the M2 wave with 65400.07 (m/strain), while O1 wave with 75700.15 (m/strain).The basic information of this well is given in Table1. 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.


## 2.2. Tidal Response Analyses

In this study, we used both O1 and M2 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 O1 and M2 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 M2 and O1 phases had opposite signs for the Gaocun well (Figure 3). A negative phase response occurred with the M2 wave, whereas a positive phase response occurred with the O1 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 Mw 7.9 Wenchuan earthquake or 2011 Mw 9.1 Tohoku earthquake).


## 2.3. Barometric Response Analyses

Figure4 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.Figure 4
Coherence coefficient of the barometric pressure and the processed water level. Data series of January 1st, 2017, to December 31th, 2017, were used to do coherence analysis.The barometric response of groundwater is obtained by cross-spectrum analysis of the water-level and the barometric-pressure time series. The time series of water-level was preprocessed to remove the effects of precipitation, tides, and recharge before the analysis [16].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].


## 3. 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:
(1)AM2′AO1′=absiωM2S/iωM2S+K′/b′ξabsiωO1S/iωO1Saqu+K′/b′ξ,ηM2=argiωM2SiωM2S+K′/b′ξ,ηO1=argiωO1SiωO1S+K′/b′ξ. Here,
(2)ξ=1+rcrw2iωrw2TβK0βrwK1βrw,β=K′Tb′+iωST1/2.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, rc is the case radius, and rw is the well radius. The subscripts M2 and O1 indicate the issues corresponding to the M2/O1 tidal frequencies.Hence, 3 input parameters (AM2′/AO1′, η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 lithology 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 Figure5. Finally, we obtained three parameters from 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].Figure 5
The analytical solution of the tidal response leaky model [9] for well Gaocun. According to the M2 and O1 phase shifts of this well, the most related analytical solution figure is chosen. ∗T, S, and K′ are in accordance with the explanations in the succeeding figure.Figure 6
Inverted vertical hydraulic conductivityK′ of the aquitard and the horizontal transmissivity T and storativity S of the aquifer from the leaky model of Wang et al. [9] for well Gaocun. The vertical dashed lines represent the start time of the 2008 Wenchuan and the 2011 Tohoku huge earthquakes, which caused no impact on those parameters. Some kick points induced by the unstable fitting perhaps incurred by the simplifications of the model assumptions (see Section 5.2) are neglected and not shown in the figure, which could disturb and obscure the dots in the figure. ∗We calculated the errors based on the data from January 1981 to December 1990 (which is relatively stable and continuous) and obtained the error of the hydraulic conductivity K′ with 1.20×10−7 (m/s), that of the horizontal transmissivity T with 7.75×10−7 (m2/s) and the horizontal storativity S with 9.31×10−6.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.


### 3.2. 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).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.
(a)
(b)Before the calculation of the barometric response, the storage coefficient of aquiferSaqu and the storage coefficient of aquitard Scon need to be set in advance. The same as the dealing of Sun and Xiang [17], we also assigned Scon=10−4 [13, 15] and Saqu=nγβLaqu/BE [17, 28] calculated in 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×103N/m3), β is the compressibility of water (4.8×10−10m2/N), and Laqu 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).Table 2
Parameter preset for the fitting of the barometric response with model C.
Model
Storage coefficient of aquiferSaqu
Storage coefficient of aquitardScon
Porosity of aquifer,n
BE
C
1.21E-03
1.00E-04
0.20
0.60


## 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:
(1)AM2′AO1′=absiωM2S/iωM2S+K′/b′ξabsiωO1S/iωO1Saqu+K′/b′ξ,ηM2=argiωM2SiωM2S+K′/b′ξ,ηO1=argiωO1SiωO1S+K′/b′ξ. Here,
(2)ξ=1+rcrw2iωrw2TβK0βrwK1βrw,β=K′Tb′+iωST1/2.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, rc is the case radius, and rw is the well radius. The subscripts M2 and O1 indicate the issues corresponding to the M2/O1 tidal frequencies.Hence, 3 input parameters (AM2′/AO1′, η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 lithology 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 Figure5. Finally, we obtained three parameters from 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].Figure 5
The analytical solution of the tidal response leaky model [9] for well Gaocun. According to the M2 and O1 phase shifts of this well, the most related analytical solution figure is chosen. ∗T, S, and K′ are in accordance with the explanations in the succeeding figure.Figure 6
Inverted vertical hydraulic conductivityK′ of the aquitard and the horizontal transmissivity T and storativity S of the aquifer from the leaky model of Wang et al. [9] for well Gaocun. The vertical dashed lines represent the start time of the 2008 Wenchuan and the 2011 Tohoku huge earthquakes, which caused no impact on those parameters. Some kick points induced by the unstable fitting perhaps incurred by the simplifications of the model assumptions (see Section 5.2) are neglected and not shown in the figure, which could disturb and obscure the dots in the figure. ∗We calculated the errors based on the data from January 1981 to December 1990 (which is relatively stable and continuous) and obtained the error of the hydraulic conductivity K′ with 1.20×10−7 (m/s), that of the horizontal transmissivity T with 7.75×10−7 (m2/s) and the horizontal storativity S with 9.31×10−6.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.


## 3.2. 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).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.
(a)
(b)Before the calculation of the barometric response, the storage coefficient of aquiferSaqu and the storage coefficient of aquitard Scon need to be set in advance. The same as the dealing of Sun and Xiang [17], we also assigned Scon=10−4 [13, 15] and Saqu=nγβLaqu/BE [17, 28] calculated in 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×103N/m3), β is the compressibility of water (4.8×10−10m2/N), and Laqu 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).Table 2
Parameter preset for the fitting of the barometric response with model C.
Model
Storage coefficient of aquiferSaqu
Storage coefficient of aquitardScon
Porosity of aquifer,n
BE
C
1.21E-03
1.00E-04
0.20
0.60


## 4. Opposite Signs for the Phase Shifts of the M2 and O1 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 M2 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 M2 and O1 responses (Figures 3 and 5) indicate substantial leakage (K′: 10-7~10-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 O1 and M2 tides also display phase shifts with opposite signs. Similarly, the opposite signs for the phase shifts of the M2 and O1 tides can be explained by the fracture model of Bower [32], where the M2 and O1 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.


## 5. Comparison between the Tidal/Barometric Response Leaky Models

### 5.1. Ideal Similar Assumptions of Those 2 Models

The tidal/barometric response leaky models are the 2-dimension model with the vertical flow in the confining layer and the horizontal flow in the aquifer (Figure8). 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).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.However, in fact, some wells could not strictly meet with the rigorous assumptions and only could do rough estimations.


### 5.2. 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.


### 5.3. 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 conductivityK′ of the confining layer and the horizontal transmissivity T of the aquifer (Tables 3 and 4).Table 3
Derived parameters from the mixed flow model.
Time period (year)
Vertical hydraulic conductivity of aquitardK′ (m/s)
Horizontal transmissivity of aquiferT (m2/s)
Storage coefficient of aquiferSaqu
Average
Standard deviation
Average
Standard deviation
Average
Standard deviation
1981-1991
4.76E-07
1.14E-07
4.73E-06
4.52E-07
7.32E-06
2.29E-06
2011-2021
4.35E-07
1.11E-07
6.97E-06
4.97E-07
5.75E-06
1.35E-06Table 4
Derived key parameters from model C of the barometric response.
Model
Horizontal transmissivity of aquiferT (m2/s)
Vertical hydraulic conductivityK′ (m/s)
Q
/
W (see [15])
C
1.11E-06
6.38E-08
354.81For groundwater response to tidal strain, we mainly used the diurnal (O1) (~0.93 cpd) and semidiurnal (M2) (~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.


### 5.4. Merits and Disadvantages of Those 2 Models

From the leaky model with tidal response, a sequence ofS, 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 aquiferSaqu and aquitard Scon 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 Saqu 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] (Saqu~10-6 (Table 3)) is much smaller than the laboratory obtained value Saqu=10−4 to 10-3 (Ss=10−6 for limestone and Laqu=717.3m 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 Saqu and Scon), and time continuity (could obtain time-varied 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.


## 5.1. Ideal Similar Assumptions of Those 2 Models

The tidal/barometric response leaky models are the 2-dimension model with the vertical flow in the confining layer and the horizontal flow in the aquifer (Figure8). 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).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.However, in fact, some wells could not strictly meet with the rigorous assumptions and only could do rough estimations.


## 5.2. 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.


## 5.3. 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 conductivityK′ of the confining layer and the horizontal transmissivity T of the aquifer (Tables 3 and 4).Table 3
Derived parameters from the mixed flow model.
Time period (year)
Vertical hydraulic conductivity of aquitardK′ (m/s)
Horizontal transmissivity of aquiferT (m2/s)
Storage coefficient of aquiferSaqu
Average
Standard deviation
Average
Standard deviation
Average
Standard deviation
1981-1991
4.76E-07
1.14E-07
4.73E-06
4.52E-07
7.32E-06
2.29E-06
2011-2021
4.35E-07
1.11E-07
6.97E-06
4.97E-07
5.75E-06
1.35E-06Table 4
Derived key parameters from model C of the barometric response.
Model
Horizontal transmissivity of aquiferT (m2/s)
Vertical hydraulic conductivityK′ (m/s)
Q
/
W (see [15])
C
1.11E-06
6.38E-08
354.81For groundwater response to tidal strain, we mainly used the diurnal (O1) (~0.93 cpd) and semidiurnal (M2) (~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.


## 5.4. Merits and Disadvantages of Those 2 Models

From the leaky model with tidal response, a sequence ofS, 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 aquiferSaqu and aquitard Scon 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 Saqu 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] (Saqu~10-6 (Table 3)) is much smaller than the laboratory obtained value Saqu=10−4 to 10-3 (Ss=10−6 for limestone and Laqu=717.3m 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 Saqu and Scon), and time continuity (could obtain time-varied 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.


## 6. 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 Figure9, 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.Figure 9
Precipitation amount of Tianjin City. Data comes fromTianjin Water Resources Bulletin (Ministry of Water Resources website), and only since the year 2000, there are available records online.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.Figure 10
Subsidence recorded by 4 GPS stations (JIXN, TJBD, TJBH, and TJWQ) within one GPS station (HECX) adjacent to the Tianjin City from the year 2011 to 2020 [37].


### 6.2. 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 O1 (diurnal) and the M2 (semidiurnal) tidal waves, through which the information of groundwater exploitation has been deleted already, except for those incidents with the same frequencies as the M2 and O1 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.


### 6.3. 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). Figure11 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.Figure 11
The original water level of well GC and well BD (the above shallower layer) during the years 2002~2015 shows a decreasing pattern. They-coordinate numbers indicate the distance between the water head and the ground.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 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.


### 6.4. Cautions for the Application of Those 2 Models (Tidal Response Leaky Model and Barometric Response Leaky Model)

Caution also should be exercised in the interpretation of the inverted hydraulic parameters with both of the 2 models because both have set numerous ideal simplified assumptions (Ideal Similar Assumptions of Those 2 Models). Those ideal simplifications might induce some inversions to be unstable, such as kick points in Figure6 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.


## 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 Figure9, 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.Figure 9
Precipitation amount of Tianjin City. Data comes fromTianjin Water Resources Bulletin (Ministry of Water Resources website), and only since the year 2000, there are available records online.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.Figure 10
Subsidence recorded by 4 GPS stations (JIXN, TJBD, TJBH, and TJWQ) within one GPS station (HECX) adjacent to the Tianjin City from the year 2011 to 2020 [37].


## 6.2. 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 O1 (diurnal) and the M2 (semidiurnal) tidal waves, through which the information of groundwater exploitation has been deleted already, except for those incidents with the same frequencies as the M2 and O1 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.


## 6.3. 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). Figure11 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.Figure 11
The original water level of well GC and well BD (the above shallower layer) during the years 2002~2015 shows a decreasing pattern. They-coordinate numbers indicate the distance between the water head and the ground.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 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.


## 6.4. Cautions for the Application of Those 2 Models (Tidal Response Leaky Model and Barometric Response Leaky Model)

Caution also should be exercised in the interpretation of the inverted hydraulic parameters with both of the 2 models because both have set numerous ideal simplified assumptions (Ideal Similar Assumptions of Those 2 Models). Those ideal simplifications might induce some inversions to be unstable, such as kick points in Figure6 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.


## 7. 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 obtainedK′ 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.



---

*Source: 4419310-2022-11-16.xml*