The origin of hydrological responses following earthquakes in a confined aquifer: insight from water level, flow rate, and temperature observations

Although many mechanisms of earthquake-induced hydrological response have been proposed in recent decades, the origins of these responses remain enigmatic, and a quantitative understanding of them is lacking. In this study, we quantitatively analyze the mechanism of coseismic response in water level and flow rate from an artesian well in southwestern China before and after multiple earthquakes and reveal the origin of the earthquake-induced hydrological response based on the monitoring data of water temperature. Water level and temperature always show coseismic step-like increases following earthquakes, which are independent of the earthquakes' epicentral distances and magnitudes. Tidal analysis finds changes in aquifer and aquitard permeability following these earthquakes, which corresponds to the post-seismic total discharge of 85–273 m³ in 20 d after earthquakes. Furthermore, we couple the flow rate and temperature data to model the mixing processes that occurred following each earthquake. The results indicate that coseismic temperature changes are the result of the mixing of different volumes of water from shallow and deep aquifers, with the mixing ratio varying according to each earthquake.

1 Introduction

Earthquake-induced hydrological anomalous phenomena, such as changes in water levels (Hosono et al., 2019; Weaver et al., 2019; H. Zhang et al., 2019, 2021), spring discharge (Manga and Rowland, 2009; Manga et al., 2016; Yan et al., 2020; Mohr et al., 2017), water temperature (Miyakoshi et al., 2020; Chien et al., 2020; Wang et al., 2012), and chemical composition (Shi et al., 2020; Skelton et al., 2014a; Kim et al., 2019; Igarashi et al., 1995; Skelton et al., 2019; Li et al., 2019; Zhang et al., 2020a) have been widely documented. Mechanisms proposed to explain these phenomena include (1) water expulsion by coseismic static strain, (2) water release from shallow crustal areas due to increased permeability following earthquakes, and (3) increased water recharge from a storage source due to an increased hydraulic head gradient. However, most of these mechanisms have been proposed depending either on observations of coseismic streamflow or water-level responses, which makes it difficult to identify the origins of the release of coseismic groundwater. These mechanisms usually do not quantitatively explain the magnitude of coseismic responses. Furthermore, the absence of integrated multi-parameter explanations for the proposed mechanisms may lead to inconsistent or even contradictory explanations of different coseismic hydrological phenomena. Although several studies have documented earthquake-induced changes in groundwater level, flow rate, and temperature, quantitative couplings of earthquake-induced groundwater flow and temperature are lacking, resulting in ambiguity in identification of the origins of earthquake-induced hydrological responses (Furuya and Shimamura, 1988; Anderson, 2010; Wang and Manga, 2010).

Based on the coseismic response of water levels, flow rate, and water temperature responses to several earthquakes in an artesian well (Dazhai well), this study aims to (1) quantitatively evaluate the mechanism of anomalous hydrological changes and (2) reveal the origin of the earthquake-induced hydrological response based on water temperature monitoring data.

2 Geological setting and data sources

2.1 Regional geological setting

Dazhai well (DZ well) is located in the Yixiangba sub-basin, 11 km southeast of the city of Pu'er in western Yunnan Province, Southwest China. The Yixiangba sub-basin belongs to the Lanping–Simao basin. The climatic characteristics in the Simao basin constitute moderate climate with four distinct seasons. The annual precipitation is about 1500 mm, but most precipitation occurs between April and October. The Simao basin is a rift basin situated in the southern section of the Sanjiang fold system. It is located in the southern Yunnan–Tibet geothermal zone, one of the most active areas of tectonic activity in China. The Lanping–Simao basin is underlain by strata ranging in age from Proterozoic to Cenozoic. The outcropping strata are related to the division of tectonic units of the basin. Due to the later denudation, the Upper and Middle Paleozoic limestones outcrop in the central and marginal uplifts of the basin. The Proterozoic metamorphic rocks occur along the Jinshajiang–Ailaoshan fault in the east part of the basin.

Tectonically, the Lancangjiang fault and Jinshajiang–Ailaoshan deep faults are the western and eastern boundaries of the Lanping–Simao basin (Deng et al., 2020). The Wuliangshan fault developed in the central basin, which is composed of several secondary faults, including Mohei, Pu'er East, Pu'er West, Zhenyuan–Puwen East, and Zhenyuan–Puwen West faults. These secondary faults are divided into several segments by a series of NE-trending faults. Fault zones can modulate fluid circulation within the shallow crust acting as conduits for groundwater flow. Veins and microfractures along faults caused by the crust dilation and contraction allow for fluid redistribution and mixing from different reservoirs (Smeraglia et al., 2018). A series of faults control the tectonic evolution of the rifting basin and also provide favorable conditions for groundwater storage and flow (Y. Zhang et al., 2019). As is shown in Fig. 1, the aquifers are classified into five categories, including loose rock aquifers, intrusive rock aquifers, metamorphic rock, clastic rock aquifers, and carbonate rock aquifers. Most of the aquifer types in the basin belong to clastic rock aquifers. The loose rock and carbonate rock aquifers occur sporadically in the northern and southern part of the basin respectively. The metamorphic rock and intrusive rock aquifers are widely distributed in the east part of the basin.

Figure 1 The hydrogeological setting of the Lanping–Simao basin. The red star indicates the distribution of 14 earthquakes, the black triangle indicates DZ well, and the solid red lines indicate faults. F₄: the Wuliangshan fault. F${}_{\mathrm{4}}^{\mathrm{1}}$: the Zhenyuan–Puwen West fault. F${}_{\mathrm{4}}^{\mathrm{2}}$: the Zhenyuan–Puwen East fault. F${}_{\mathrm{4}}^{\mathrm{3}}$: the Pu'er West fault. F${}_{\mathrm{4}}^{\mathrm{4}}$: the Pu'er East fault. F${}_{\mathrm{4}}^{\mathrm{5}}$: the Mohei fault.

2.2 Information of monitoring well

DZ well (constructed in 1984 and upgraded to a digitized observation system in 2001) is located in fractured sandstone–siltstone rock in the Jinmenkou–Shigaojing part of Pu'er West fault (Fig. 1). DZ well is an artesian well with a discharge rate of 0.3–0.6  L s⁻¹, a depth of 112.27 m with a casing diameter of 150 mm, a blank casing from 0–26.69 m with an open section between 26.69–112.27 m, and diameters of 120 and 110 mm at depths of 26.69–83.7 and 83.7–112.27 m, respectively. The elevation level of Dazhai well is about 1471 m. Groundwater level, flow rate, and water temperature are observed in DZ well. The lithology of the wellbore in DZ well is as follows (Fig. 2): 0–7.51 m is quaternary deposit, 7.51–14.03 m is argillaceous siltstone with strong weathering, 14.03–26.69 m is silty mudstone, 26.69–78.90 m is argillaceous siltstone, 78.90–110.76 m is detrital sandstone with fissure developed, and 110.76–112.27 m is mudstone. In order to understand the deeper lithology and the distribution of aquifers, we also collected the aquifer lithology and wellbore structure of Dazhai deep well (Fig. S1 in the Supplement), which is drilled with a depth of 410 m and only 50 m away from Dazhai well. Both wells are artesian wells. Sandstone has the characteristics of high mechanical strength and permeability and is considered the most important type of aquifer in this area. Conversely, the characteristics of mudstone are low mechanical strength and permeability, and it is considered an aquitard. In combining these two wellbores, we found that the sandstone stratum is separated by mudstone at depth between 105–117 m, forming two confined aquifers, and the groundwater head in the underlying sandstone aquifer (1480 m) is higher than the shallow sandstone aquifer (1475 m). In addition, the hydrochemical types of Dazhai well and Dazhai deep well are HCO₃–Na ⋅ Ca and HCO₃–Ca ⋅ Na, respectively. There is no significant difference in groundwater hydrochemical type in samples from the shallow and deep aquifers.

Figure 2 Aquifer lithology and specifications of the DZ well aquifer system.

2.3 Data collection

Earthquake groundwater monitoring wells in China are categorized as two types: (1) static water-level observation for the non-artesian well and (2) dynamic water-level observation for the artesian well (CEA, 2012). DZ well is an artesian aquifer, and the hydraulic head is always above the surface. Thus, dynamic water-level observation is employed. Such a groundwater monitoring system has been widely used in Chinese earthquake groundwater monitoring wells in artesian aquifers. Dynamic water-level observation wells tend to have better confinement than the non-artesian flowing wells. Previous studies revealed that three external influence factors affect the sensitivity of the earthquake-induced hydrogeological responses, including peak ground velocity (PGV), aquifer confinement, and well location relative to local faults (Zhang et al., 2021). In additional, many studies have proposed that the tidal signals from aquifers with different confinement are different. The aquifers with better confinement could record clear and large magnitude of M₂ tidal wave, while the aquifers with less confinement may show weak tidal signals (Rahi, 2010; Turnadge et al., 2019; Zhang et al., 2021). Thus, dynamic water-level observations tend to show more sensitivity to earthquakes. In order to observe the variation of water-level fluctuation, the wellbore is constructed above the ground and equipped with the discharge port (Fig. 2). The purpose of setting the discharge port is to discharge a part of groundwater and keep the water-level fluctuation in a certain range. Meanwhile, the discharge can be measured through this port sing the stopwatch capacity method; i.e., the required time per unit volume of water is measured. Without the discharge port, the hydraulic head would be 4 m above the surface, which is difficult for monitoring. A piezometric tube is used to calibrate the data of water level. For such a monitoring system, the water level refers to the distance between the water table and the center line of the discharge port. Dynamic water-level observations in DZ well are recorded hourly with a resolution of 1 mm using a LN-3A pressure transducer (developed by the Institute of Earthquake Forecasting, China Earthquake Administration, Beijing, China). Water temperatures are digitally recorded every minute using a quartz crystal thermometer (SZW-1A, developed by the Institute of Crustal Dynamics, China Earthquake Administration, Beijing, China) with a resolution of 0.1 mK (Shi and Wang, 2014). The temperature probe is set at a depth of 100 m and removed from the influences of other factors such as air temperature, precipitation, and pumping. Monitoring data are transferred to and stored at the Yunnan provincial earthquake administration, Kunming, China.

We collected water-level, flow rate, and water temperature data for the period 2004–2015. Water levels show seasonal fluctuations with a higher level during the summer wet season and a lower level during the drier winter season. The δD and δ¹⁸O values of groundwater from Dazhai well are −80.63 ‰ and −11.17 ‰, respectively, which are close to the global meteoric water line (GMWL): δD = 8δ¹⁸O + 10 (Yang et al., 2005). The stable-isotope data indicate that the recharge of DZ well aquifer is primarily from meteoric waters. According to the relationship between recharge elevation and δD in the eastern Tibetan Plateau δD = −0.026 H (m) − 30.2 (Yu, 1997), the recharge elevation of groundwater from Dazhai well is about 1940 m. The estimated recharge elevation is close to the elevation of mountains nearby Dazhai well.

3 Observation of coseismic water-level and temperature changes

Three patterns of water-level coseismic response have been identified based on the distance between the length of the ruptured fault and the epicentral distance: (1) step-like changes in near-field earthquakes (< 1 ruptured fault length); (2) gradual changes, both increases and decreases, at intermediate distance (1–10 ruptured fault length); and (3) oscillations in far-field earthquakes (> 10 ruptured fault length). During the 2004 to 2015 monitoring intervals, these responses were observed for 14 earthquakes, including the 2004 M_w 9.0 Sumatra earthquake, 2008 M_w 7.9 Wenchuan earthquake, and 2015 M_w 8.2 Nepal earthquake (Fig. 1). There is no barometric pressure monitoring data of 2004–2015 in the monitoring site of DZ well; thus no barometric pressure correction had been done for the water level. The raw water-level data are shown in Fig. S2. We have removed the effect of Earth tide in calculating the coseismic water-level changes. The air temperature has no effect on the variation of water (Sect. S3 in the Supplement). The coseismic and post-seismic changes in water level and temperature can be clearly identified (Fig. 3) and show sharp step-like rises, with anomalous changes following earthquakes lasting approximately a month (summarized in Table 1).

Table 1 Seismic properties for 14 earthquakes and their effects on water levels and temperatures in the Dazhai well, western Yunnan Province, China.

Figure 3 Water-level and temperature changes induced by 14 earthquakes in Dazhai well, western Yunnan Province, China. The water-level changes have been corrected for Earth tide. Vertical dashed lines indicate earthquake events.

All water-level coseismic responses show step-like increases, with magnitudes ranging from a minimum of 58 mm to a maximum of 410 mm caused by the 3 June 2007 M_w 6.4 (EQ4). Water-level coseismic changes reached 228 mm in response to epicentral distances greater than 3000 km (EQ5). Water temperature induced by earthquakes ranged from 2.1 to 24.5 mK, where the maximum was caused by the 24 March 2011 M_w 7.2 earthquake (EQ7).

Water levels also show the influence of Earth tides, with a maximum variation induced by semidiurnal and diurnal tides being < 9 mm (Yang et al., 2005), which is less than the step-like changes induced by earthquakes. Seasonal fluctuations in water level from rainy to dry seasons were < 200 mm and showed no long-term trend. All earthquake-induced groundwater-level responses were greater than 200 mm, which cannot be attributed to Earth solid tides and meteorological factors.

Wang (2007) introduced the concept of seismic energy density to quantify the magnitude and epicentral distance of earthquakes to build the relationship between seismic activity and hydrologic response, where a threshold value of 10⁻³ J m⁻³ was concluded to be sufficient to trigger a hydrological coseismic response. The seismic energy density can be calculated by the following equation:

$$\begin{array}{}\text{(1)}& \log r=\mathrm{0.48}{M}_{\mathrm{w}}-\mathrm{0.33}\log e-\mathrm{1.4},\end{array}$$

where r is the epicentral distance, M_w is the magnitude of the earthquake, and e (J m⁻³) is the seismic energy density. Table 1 summarizes the seismic energy density of the 14 earthquakes, which ranged from 0.005 to 2.242 J m⁻³. Detailed information is summarized in Table 1. Those earthquakes that showed hydrological response and no hydrological response are plotted in Fig. S4. As is shown in Fig. S4, although the seismic energy density of most earthquakes is < 10⁻³ J m⁻³, the seismic energy density of a few earthquakes ranging from 10⁻³ to 10⁻² J m⁻³ causes no hydrological response. The minimum seismic energy density required to trigger water-level and temperature coseismic changes was 0.0054 and 0.005 J m⁻³, respectively, where both values are > 10⁻³ J m⁻³, approaching the empirical threshold for triggering hydraulic property changes, which means that this well is highly sensitive to seismic activity. We next explore the mechanism of these coseismic responses and their quantitative assessment and then identify the origin of excess recharge to the aquifer induced by earthquakes.

4 Results

4.1 Changes in horizontal hydraulic diffusivity

Aquifer parameters such as transmissivity, storage coefficient, and hydraulic diffusivity play important roles in groundwater flow. During seismic activity, anomalous changes in water level may be attributed to earthquake-induced changes in aquifer parameters. The tidal analysis method establishes a direct connection between hydrological responses and aquifer parameter changes, which is also an effective method for the long-term monitoring of variation in aquifer parameters (Hsieh et al., 1987; Xue et al., 2013; Shi and Wang, 2016). Water-level changes induced by Earth solid tides in a wellbore show semi-diurnal and diurnal fluctuations, which are controlled by aquifer properties and the wellbore structure. The amplitude is the ratio of the water-level amplitude to the theoretical volumetric strain, and the phase shift is the time lag between the Earth solid tide fluctuation and water-level variation, which differ because of the time required for groundwater to equilibrate the well. Both can be used to calculate the transmissivity and storage coefficient:

$$\begin{array}{}\text{(2)}& {\displaystyle }X/H=(E^{\mathrm{2}}+F^{\mathrm{2}}{)}^{-\mathrm{1}/\mathrm{2}}\text{(3)}& {\displaystyle }\mathit{\eta }=-{\tan }^{-\mathrm{1}}(F/E)\text{(4)}& {\displaystyle }E\approx \mathrm{1}-{\displaystyle \frac{\mathit{\omega }{r}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}T}}\mathrm{Kei}\left(\mathit{\alpha }\right),F\approx {\displaystyle \frac{\mathit{\omega }{r}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}T}}\mathrm{Ker}\left(\mathit{\alpha }\right),{\mathit{\alpha }}_w=({\displaystyle \frac{\mathit{\omega }S}{T}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}{r}_{\mathrm{w}},\end{array}$$

where $X/H$is the amplitude; X is the amplitude of tide of the water level in the well; H is the amplitude of the fluctuating pressure head in the elastic aquifer responding to the tidal stress; η is the phase shift; T is transmissivity (m² d⁻¹); S is the storativity (dimensionless); Ker and Kei are zero-order Kelvin functions; r_w and r_c are the radii of the well and well casing, respectively; and ω is the tidal frequency (Hz).

Figure 4 Phase shift and tidal factor of M₂ tidal constituent using water-level observations in Dazhai well between 2004 to 2015. Vertical dashed lines indicate earthquake events.

Spectral analysis was conducted to identify the multiple components inducing water-level fluctuations in DZ well (Fig. S5). The obvious amplitudes were observed in a frequency band of 0–3 cpd, indicating that the water level in Dazhai well responds to Earth tides (Qu et al., 2020). In order to obtain the phase shift and amplitude information of the water-level time series, we used tidal analysis in the calculations. We used a common tidal analysis program (Baytap-G) and set a period of 30 d with an overlap of 15 d (Tamura et al., 1991). For Earth tide, O₁, K₁, M₂, and S₂ are the main waves that account for at least 95 % of all tidal waves. Here, we only chose the M₂ wave for further analysis as it is stable and relatively unaffected by barometric pressure. We focus on the evolution in aquifer properties induced by seismic activity based on post-seismic water level to recovery to the pre-seismic value that ranges from 13 to 70 d. We considered earthquakes that occurred within 3 months of each other as a single earthquake to reduce post-seismic effects on the aquifer properties. During this time interval, earthquakes with the greatest seismic energy density were selected for tidal analysis, which totaled seven earthquakes (EQ1, EQ4, EQ6–EQ8, and EQ11–12). The 3-month averages of phase shift and amplitude before and after earthquakes were used to calculate the aquifer parameters. We collect the monitoring data of barometric pressure in the city of Simao and find that barometric pressure fluctuation would not have an effect on the result of our tidal analysis (Figs. S6, S7 and S8). The phase shift and tidal factor changes following earthquakes can be clearly identified in Fig. 4, which are summarized in Table 2.

Table 2 Aquifer parameters estimated using tidal analyses before and after selected earthquakes in Dazhai well, western Yunnan, China.

As is shown in Fig. 4, the errors are much smaller than the changes that are caused by the earthquakes, which represent the 95 % confidence interval. There is no periodic variation on the tidal factor. During the monitoring period, the tidal factors range from $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{3}}$, and the average value is $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$. During the aseismic period, the rate of changes ranges from 0 to 75% compared with the mean value of tidal factor. However, the changes are 50 %, 237.5 %, 87.5 %, 150 %, 12.5 %, 100 %, and 25 % following EQ1, EQ4, EQ6, EQ7, EQ8, EQ11, and EQ12, which are greater than that during the aseismic period. Thus, the changes of tidal factor induced by earthquakes are real and significant.

The average value of phase shift is −20^∘ during 2004–2015, and it ranges from −31 to −11^∘ during the period without earthquakes. However, the phase shifts following EQ1, EQ4, EQ7, and EQ8 are greater than −10^∘, which is also significant with that during aseismic period. During the pre-seismic and post-seismic period of EQ11 and EQ12, the phase shift following these earthquakes is also > −10^∘, and it shows a wide range of fluctuation which looks like “noise”, but it is not actually noise. The reason for the phenomenon may be that the epicentral distance of EQ11 and EQ12 is 4 and 100 km, respectively, and the aquifer properties are continuous change caused by crustal strain change during the earthquake preparation process. In addition, because of the short interval between the two consecutive earthquakes, the aquifer properties may have been affected by the previous earthquake and not recovered.

As proposed by Hsieh et al. (1987), the phase shift and tidal amplitude are related to the aquifer transmissivity and storage coefficient, so that changes in phase shift and tidal amplitude also indicate changes in aquifer hydraulic properties following earthquakes. The tidal amplitude reflects the ability of an aquifer system to respond to Earth tidal changes induced by earthquakes and is affected by aquifer properties (e.g., undrained compressibility, Poisson's ratio, and fluid compressibility), while variations in phase shift are related to permeability changes (i.e., small phase shifts are related to high aquifer permeability, while large ones are related to low permeability). The hydraulic diffusivity D (m² s⁻¹; equal to the ratio of T to S), is introduced to accurately depict the changes in aquifer parameters, where results imply that earthquakes cause the change of hydraulic diffusivity ranges from 9 % to 515 %.

4.2 Changes in vertical permeability

Hsieh's model only analyzes the earthquake-induced horizontal permeability changes. The vertical permeability change caused by earthquakes can be detected by comparing the amplitude and the tidal phase shift (Liao et al., 2016; Wang et al., 2016). As is shown in Fig. 5, the amplitude and phase shift inferred from before and after earthquakes fall into two distinct groups (red dots and blue dots). In addition, for a given value of H, the X versus η relationship can be represented by a single curve in the plot, which is uniquely determined by T and S. Dots before earthquakes (blue) and after earthquakes (red) fit two curves with different values of H, which indicated the changes of vertical permeability in the aquifer system (Liao et al., 2016).

Figure 5 Amplitude versus phase shift at the frequency of the M₂ wave. Dots before earthquakes (blue dots) and after earthquakes (red dots) fit two curves that represent the theoretical model (Eqs. 2 and 3) with different values of H, respectively.

5 Discussion

5.1 Mechanisms of water-level coseismic responses

Static and dynamic strain mechanisms have been proposed to explain coseismic water-level responses where static strain changes aquifer permeabilities, resulting in variations in pore pressure and fracture aperture (flow chart of different mechanisms and methods in Fig. S9), while the effects of dynamic strain include (1) undrained consolidation and liquefaction, (2) drainage from the unsaturated zone, (3) bubbles mobilization through highly fracture, (4) the formation of new fractures, and (5) unclogging of preexisting fractures due to the propagation of seismic waves. In this section, we will discuss the impacts of static and dynamic strain on water-level coseismic responses in DZ well.

5.1.1 Static strain

Static strain induced by fault slip may lead to the dilation of fractures which will change the pore pressure and the jog volume. For our data, static strain produced by large, distant earthquakes with epicentral distances > 2000 km (EQ1 and EQ8) is lower than for earthquakes listed in Table 2; thus we use the near-field earthquakes listed in Table 2 to analyze the change of jog volume and fluid pressure.

The elastic or permanent deformation of crust associated with the fault slip causes static strain. It is assumed that the fault jog was occupied by the fluid before earthquakes. When the frictional strength of the fault is exceeded under the accumulation of tectonic stress, we may estimate the relative volume change of fault jog according the seismic moment (M₀),

$$\begin{array}{}\text{(5)}& M_{\mathrm{0}}=\mathit{\mu }As,\end{array}$$

where μ is the shear rigidity of the host rock, A is the surface area of the rupture, and s is the average slip over the rupture area. The seismic moment can be converted to the moment magnitude (M_w):

$$\begin{array}{}\text{(6)}& M_{\mathrm{w}}=\mathrm{2}/\mathrm{3}\log M_{\mathrm{0}}-\mathrm{10.7}.\end{array}$$

The rupture length (L) of seismogenic fault is the function of moment magnitude (Wells and Coppersmith, 1994), as follows:

$$\begin{array}{}\text{(7)}& \log L=\mathrm{0.74}{M}_{\mathrm{w}}-\mathrm{3.55}.\end{array}$$

The down-dip width (W) and the average slip (s) can be calculated by the established empirical relation (Weatherley and Henley, 2013):

$$\begin{array}{}\text{(8)}& {\displaystyle }W=L{M}_{\mathrm{w}}<\mathrm{4.0}\text{(9)}& {\displaystyle }W=\mathrm{15.0}{L}^{\mathrm{2}/\mathrm{3}}\mathrm{4.0}<M_{\mathrm{w}}<\mathrm{6.0}\text{(10)}& {\displaystyle }W=\mathrm{17}\mathrm{km}{M}_{\mathrm{w}}>\mathrm{6.0}\text{(11)}& {\displaystyle }\log s=\mathrm{0.833}\log L-\mathrm{3.84}.\end{array}$$

Based on the above parameters, the relative volume change of fault jog could be estimated by

$$\begin{array}{}\text{(12)}& \mathrm{\Delta }V=w{D}_{j}s,\end{array}$$

where w is the vertical extent of the jog, and D_j is the jog length; it is assumed that D_j= 1 m (Weatherley and Henley, 2013). The relative volume change of fault jog is summarized in Table 3 following different near-field earthquakes.

Table 3 The relative change of volume and fluid pressure following the near-field earthquakes.

Based on the assumption of isothermal expansion, the fluid pressure changes in the fault jog associated with the volume change of fault jog are considered the behavior of ideal gas with the equation of state PV = nRT (Weatherley and Henley, 2013). According to the 11-year observed mean value of water level, the initial pressure (P) of the fluid before the earthquake is calculated as

$$\begin{array}{}\text{(13)}& P\left(z\right)=\mathit{\rho }gz,\end{array}$$

where ρ= 1000 kg m⁻³, g= 9.81 m s⁻² is the gravitational acceleration, and z is the depth. Then, the relative change of fluid pressure is calculated by

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{P_{\mathrm{f}}}{P_{\mathrm{i}}}}={\displaystyle \frac{V_{\mathrm{i}}}{V_{\mathrm{f}}}},\end{array}$$

where the subscripts “i” and “f” indicate values before (initial) and after (final) the earthquake. The significant differences of final fluid pressure varying with the different depths are summarized in Table 3. For different earthquakes, the magnitudes of relative change of fluid pressure are different. According to P=ρgh, the relative change of fluid pressure can cause the magnitude of coseismic response of water level to exceed 80 m in theory, which is inconsistent with the observed value of DZ well. Due to the action of the discharge port on the wellbore, a part of coseismic increased groundwater is discharged by the discharge port. The actual maximum change of water-level coseismic response should be greater than the observed value. We estimate that the actual maximum change of water-level coseismic response should be 0.35–0.83 m, which is the sum of the observed maximum change and the conversion of the volume of discharge groundwater into the height of water level in wellbore, based on the equation of V=πR²H, where V is the volume of discharged groundwater, H is the height of water level, and R is the radius of wellbore. But these results are still 2 orders of magnitude smaller than the theoretical value of water-level change estimated by the relative change of fluid pressure. Thus, the water-level coseismic response may not be attributed to the fluid pressure change induced by static strain.

The relative change of fluid pressure only calculates the theoretical value of relative magnitude of water-level change but cannot judge whether water level is rising or declining following earthquakes theoretically. The Okada (1992) model is introduced to further evaluate the relationship between the water-level coseismic response and static strain, which indicates that the change of static strain causes the water levels to decline in the expansion zone and rise in the contraction zone (Okada, 1992; Toda et al., 2011). Table S1 in the Supplement summarizes fault-plane solutions obtained from the USGS for hypocentral depth, strike, dip, and rake.

Observed coseismic water-level changes following EQ12 are consistent with this model, while the coseismic responses following EQ4 and EQ7 are opposite to the prediction (Fig. S10). We conclude that the distribution of coseismic water-level responses during seismic activity is more complex than the quadrantal distribution expected from the theoretical model. Thus, based on the analysis of the relative change of fluid pressure and the Okada model, the coseismic response of water level in DZ well cannot be attributed to the static strain because it only explains a subset of the coseismic water-level responses. In addition, the large distant earthquakes (epicentral distances > 2000) also cause the coseismic responses of water level in DZ well, while static strain is extremely small for these distant earthquakes, thus implying that the mechanism of coseismic responses is attributed to dynamic strains.

5.1.2 Dynamic strain

As DZ well is located in a highly fractured aquifer, the coseismic hydrological responses are more likely caused by the mechanism of dynamic strain. Of the five potential dynamic causes, the first two (undrained consolidation and liquefaction, unsaturated zone drainage) can be omitted because the observation well is located within consolidated rocks, and no liquefaction was observed during seismic activity. Moreover, the variation of water-level in DZ well is sensitive to Earth tides, which indicates that the content of compressible gaseous phase is less than 10⁻⁴ % of the aquifer (Brodsky et al., 2003). Thus, the mechanism of bubble mobilization through high fracture is unlikely. In addition, the result of tidal analysis indicates that these earthquakes caused changes in aquifer properties; the phase shift increased following earthquakes and recovered within several weeks to months afterward, indicating that the earthquake-induced permeability changes were transient, not permanent. The change of the permeability and poroelastic properties within the fault damage zone may be attributed to the dynamic strain caused by the passage of seismic waves. Thus, we consider that earthquake-induced dynamic stress leading to unclogging of preexisting fractures under the propagation of seismic waves is a plausible mechanism for the observed behavior. According to the aquifer lithology of Dazhai deep well, there were two different aquifers separated by a thin aquitard (mudstone), and the pore pressure of the deeper aquifer is higher than the shallow one; such a hydrogeological setting makes the upwelling of groundwater from deep aquifers possible when the hydraulic connection is created following earthquakes (Fig. 6). In the unclogging process, precipitates and colloidal particles are removed from the clogged fractures by the propagation of seismic waves, which connects aquifers at different depth. After earthquakes, permeabilities gradually recover through geochemical and biogeochemical processes or by the re-clogging of fractures, causing water levels to recover to pre-earthquake levels (Manga et al., 2012; Xue et al., 2013; Shi and Wang, 2015).

Figure 6 Conceptual model for the recharge process and hydrological coseismic response in Dazhai well.

Yet, unclogging only qualitatively explains the change of water levels following earthquakes and does not explain the different magnitudes of changes following different earthquakes. Alternatively we infer that the different magnitudes of water-level change are associated with the total volume of excess recharge to the aquifer during seismic activity, where the volume of aquifer recharge can be described by Wang et al. (2004):

$$\begin{array}{}\text{(15)}& q={\displaystyle \frac{\mathrm{8}Q}{{\mathit{\pi }}^{\mathrm{2}}\mathit{\tau }}}\sum _{n=\mathrm{0}}^{\mathrm{\infty }}\exp \left[-(\mathrm{2}n+\mathrm{1}{)}^{\mathrm{2}}{\displaystyle \frac{t}{\mathit{\tau }}}\right],\end{array}$$

where Q is the total excess water recharge to the aquifer from deeper aquifer, t is the time since the earthquake occurred, and τ is characteristic time. The solution to the groundwater flow equation can be determined by two parameters (Q and τ) that are determined by fitting Eq. (15) to the observed discharge time series using the nonlinear least-squares Marquardt–Levenberg algorithm (Fig. 7). The accuracy and efficiency of model fitting results are shown in Table S2 and Fig. S11.

Figure 7 Model results compared with normalized discharge data for the discharge coseismic response in DZ well.

Excess recharge predicted by the one-dimensional conceptual model is summarized in Table 4. The excess recharge caused by earthquakes differs for different earthquakes within the range of 85–273 m³. We conclude that this is the most plausible reason for the variation in water-level change magnitudes following each earthquake.

Table 4 Parameters estimated by the data of discharge coseismic response in Dazhai well, western Yunnan Province, China.

5.2 Possible mechanism of water temperature coseismic change

The dynamic strain mechanism describes the process of coseismic water-level change, while the one-dimensional diffusion model can explain the different magnitudes of anomalous water-level change following each earthquake. Because they cannot reveal the source of water recharge to the aquifer, we analyzed the coseismic water temperature changes to further identify the origin and processes of the coseismic groundwater responses. Several mechanisms have been proposed for groundwater temperature changes following earthquakes, including heat loss through the walls of boreholes (Kitagawa and Koizumi, 2000), gas escape from groundwater (Ma, 2016; He et al., 2017), friction induced by seismic waves (Fulton et al., 2013; Tanaka et al., 2006; Mori et al., 2010), and changes in the mixing ratios of groundwaters with different temperatures (Shi et al., 2018; Mogi et al., 1989). Although heat can be propagated by radiation, convection through pores, and conduction through solid phases, solar radiation only occurs at the Earth's surface, which has little effect on groundwater temperature. Also, temperatures should change gradually if groundwater exchanges heat with the rock matrix, which is inconsistent with our observations. The gas-escape mechanism explains groundwater temperature changes caused by large amounts of dissolved gases escaping from groundwater (a phase change from liquid to gas), which absorbs heat from the surrounding environment that is then dissipated, causing the water temperature to decrease. Yet this explanation is also inconsistent with our step-like increases. Friction heat caused by seismogenic faults may explain coseismic water temperature responses in some locations close to earthquake epicenters but cannot cause temperature responses at great distances.

Because groundwater is an important carrier fluid for heat transfer, changes in water level, flow rate, and temperature occur simultaneously, which indicates that groundwater flow plays an important role in coseismic responses of water temperature for our system. Thus, the mechanism of mixing of different sources of water is likely to explain the synchronous changes in groundwater discharge and temperature. Several previous studies have documented that mixing of deeper fluids and meteoric water in the crust may play a key role in changing the physical and chemical characteristic of the groundwater (Kaown et al., 2021; Li et al., 2019; Skelton et al., 2014b; Jellicoe et al., 2022; Woith et al., 2013). The fault zone could act as a good conduit that connects the shallow and deep groundwater. For the hydrothermal spring in the fault zone, it usually mixes with the shallow cold water during its rising to the surface. And groundwater in the fault zone aquifer may also originate from different depths of the aquifers. Previous studies suggested that earthquakes could lead to the enhancement of permeability in the low permeability layer (Fischer et al., 2017) or reduction of permeability in high flow path (Shi et al., 2018), causing the changes of mixing ratio between shallow and deep groundwater, and thus lead to the changes of flow rate and chemical components in groundwater (Sato et al., 2020; Li et al., 2019). However, few studies have investigated the earthquake-induced changes both in discharge and temperature in a quantitative way. An end-member mixing model is widely used to quantify the different contributions of groundwater from two or more sources, which is based on assumptions that the spatial distribution of groundwater with distinct characteristics is nonuniform. In this study, we use the end-member mixing model to describe the anomalous change of discharge and water temperature following earthquakes. The possible end-member is analyzed and identified in the following part.

We observed an increase of the water level and flow rates which are not related with precipitation as only 3 and 0.9 mm, which occurred in December 2004 (EQ1) and February 2014 (EQ11), respectively. Wavelet coherence analysis is widely used to identify the influences on water temperature variation in the absence of earthquakes (Zhang et al., 2020a). As shown in Fig. S12 of the Supplement, there is high coherence between meteorological and hydrological data for around a year throughout the measurement period. However, the mean phase angles between discharge (water temperature) and rainfall are about 45^∘, indicating a lead of rainfall about 1–1.5 years, which implies that the excess recharge to the aquifer is not directly originated from precipitation. Although DZ well only penetrates one aquifer, according to the borehole information of Dazhai deep well, there were two different aquifers separated by a thin mudstone aquitard. Seismic activities cause the increased permeability of the aquitard and make the upwelling of groundwater from the deep aquifer possible when the hydraulic connection is created (Fig. 6). Thus, we consider the mixing process between them.

Figure 8 Observed and estimated discharge for shallow and deep source following selected earthquakes.

For the selection of end-members, there are no limitations on the number of them in the hydrologic model, but the general rule is to use as few as possible to describe the observed anomalous phenomena (Laaksoharju et al., 2008). Thus, we assume that the different magnitudes of water temperature change are also caused by the mixing of these from distinct end-members at different ratios following earthquakes. The earthquake-induced mixing process between the shallow and deep aquifer can be described by the following equations (Kitagawa and Koizumi, 2000):

$$\begin{array}{}\text{(16)}& {\displaystyle }{D}_{\mathrm{observed}}=D_{\mathrm{s}}+D_{\mathrm{d}}\text{(17)}& {\displaystyle }{T}_{\mathrm{observed}}={\displaystyle \frac{T_{\mathrm{s}}{D}_{\mathrm{s}}+T_{\mathrm{d}}{D}_{\mathrm{d}}}{D_{\mathrm{s}}+D_{\mathrm{d}}}},\end{array}$$

where D_observed is the observed discharge; D_s and D_d denote the discharge from shallow and deep aquifers, respectively; T_observed is the observed water temperature; and T_s and T_d are the water temperatures of shallow and deep aquifers, respectively. We infer that the water temperature responses caused by earthquakes are related to changes in discharge from the shallow and deep aquifers.

During the monitoring period, temperatures ranged between 21.7788 and 21.9325^∘, a difference of 153.7 mK that remained within this range even following seismic events. To explain temperature variation following these events, we employed a two-end-member mixing model, where we assume two aquifers at different depth that are hydraulically connected following earthquakes: one is shallower with a lower temperature (21.7788^∘), while the other is deeper with a higher temperature (21.9325^∘), which corresponds to the extremal temperatures observed during the period of record.

Discharge from each end-member is then estimated using

$$\begin{array}{}\text{(18)}& {\displaystyle }{D}_{\mathrm{s}}=D_{\mathrm{observed}}{\displaystyle \frac{T_{\mathrm{observed}}-T_{\mathrm{d}}}{T_{\mathrm{s}}-T_{\mathrm{d}}}}\text{(19)}& {\displaystyle }{D}_{\mathrm{d}}=D_{\mathrm{observed}}(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{observed}}-T_{\mathrm{d}}}{T_{\mathrm{s}}-T_{\mathrm{d}}}}).\end{array}$$

Figure 8 presents the calculated discharge for each end-member, with end-members 1 and 2 changing following earthquakes; for EQ1, EQ11, and EQ12, the estimated discharge for shallow aquifer remained unchanged before and after earthquakes, while discharge from deep aquifer increased sharply after earthquakes, which may be due to the permeability increase in deep aquifers, which results in more warm groundwater entering the wellbore. For the other four earthquakes, the estimated discharge from shallow and deep aquifers varied, with the magnitude of change varying according to the earthquake. The change in discharge from deep aquifers was always greater than that from shallow aquifers, which is the main cause of the coseismic water temperature response. Thus, we conclude that the coseismic water temperature response can be attributed to earthquake-induced permeability increases, which result in warmer groundwater entering the wellbore. The different magnitudes of water temperature change following each earthquake are caused by changes in the mixing ratios of shallow and deep aquifers with different temperatures.

6 Conclusion

We document groundwater levels, temperatures, and flow rate responses following 14 earthquakes in an artesian well in western Yunnan Province, China. The well always shows step-like groundwater-level responses and increased flow rates and temperatures following earthquakes. Tidal analysis indicates that the aquifer's hydraulic diffusivity changes following earthquakes. Due to the unclogging of preexisting fractures under the propagation of seismic waves, the hydrogeological setting makes the upwelling of groundwater from the deeper aquifer possible when the hydraulic connection is created following earthquakes, which is the likely mechanism explaining the coseismic water-level response. The excess groundwater recharge from the deep aquifer, ranging from 85 to 273 m³ for the 14 earthquakes, leads to different magnitudes of water-level changes. A groundwater mixing model based on observed water temperatures is used to estimate mixing processes of waters with different temperatures before and after earthquakes, which shows that earthquakes induce recharge from both shallow and deep aquifers.