Determining Stress State of Source Media with Identified Difference between Groundwater Level during Loading and Unloading Induced by Earth Tides

Abstract

The groundwater level might be adopted as a useful tool to explore pre-seismic stress change in the earth crust, because it circulates in the deep crust and should be altered by the processes associated with the preparation of earthquakes. This work makes a new attempt that applies the load/unload response ratio (LURR) technique to study the stress state of the source media before the large earthquakes by calculating the ratio between the water levels during the loading and unloading phases. The change of Coulomb failure stress induced by earth tides in the tectonically preferred slip direction on the fault surface of the mainshock is adopted for differentiating the loading and unloading periods. Using this approach, we tested the groundwater level in the wells near the epicenters of some large earthquakes that occurred in the Sichuan-Yunnan region of southwest China. Results show that the LURR time series fluctuated narrowly around 1.0 for many years and reached anomalously high peaks 2~8 months prior to the mainshocks. For the earthquakes with multiple observation wells, the magnitude of the maximum values decreases with the distance from the epicenter. The underlying physics of these changes should be caused by the pre-seismic dilatancy. The corresponding volume variations in the crust could be observed in the geodetic time series in the same neighborhoods and during the same period.

1. Introduction

The instrumentally recorded data of groundwater level changes have been widely applied to study processes associated with the occurrence of earthquakes, and reported in many studies [1,2,3]. Roeloffs [4] revealed the response of water level fluctuations in a well near the Parkfield to distant and local earthquakes. King et al. [5] detected earthquake-related changes by analyzing co-seismic water level changes at 16 closely clustered wells with different depths. Brodsky et al. [6] suggested the mechanism of groundwater level changes and pressure induced by distant earthquakes. The results seem to support the notion that observation of groundwater level changes helps us understand the responses of the hydrogeologic system to earthquakes [7].

However, the exploration of the water level fluctuations before the occurrence of earthquakes is much more difficult. The current indexes of the anomalies for detection of the time, location, and magnitude of future large earthquakes are mostly obtained by using the statistical techniques [8], e.g., the groundwater level anomalies preceding the large earthquakes in China and Japan [9], the short-term pre-seismic changes in groundwater level in Turkmenistan [10], and the precursory changes before moderate local earthquakes and some large distant earthquakes in the well water level near the Tono mine in central Japan [11]. Due to the lack of understanding of the physical mechanism of earthquake preparation, less information can be extracted from the ambiguous anomalies and their uncertain time scales for the studies of earthquake potential [12,13].

In our latest studies [14], we find that the notion of Load/unload response ratio (LURR) may be applied to explore physical processes associated with the preparation of earthquakes via groundwater level. The LURR method is proposed by Yin et al. [15] and developed by Zhang et al. [16] and Yu et al. [17]. They believe that the variation of LURR time series is a proxy for the stress state of a heterogeneous brittle system (Figure 1). When a nonlinear system is in the elastic phase, nearly the same responses (e.g., the strain, energy release, and so on) to the loading and unloading are observed, and LURR ≈ 1.0; whereas when the system is stressed beyond the elastic limit, the responses become significantly different, and LURR > 1.0. Thus, the LURR value might be served as a useful parameter to assess the criticality of earthquakes. Several publications have reported anomalous enhancements prior to a large earthquake in the time series of LURR [18,19,20].

Traditional LURR practice is based on the calculation of the ratio between Benioff strain releases during the time periods of loading and unloading, involving the evaluation of the Coulomb failure stress change [21] induced by earth tides and a regional earthquake catalog. However, recent studies have found that the tidal responses can also be observed in the well water level [22,23,24,25]. Meanwhile, the circulation of groundwater in the crust, especially in fault zones, can reach 15~20 km, corresponding to the depth of most earthquakes [26,27]. The foreshocks or cracks resulted from the pre-seismic stress accumulation can cause rock dilatancy [28], and the magnified volume in the crust may lead to changes in the groundwater level. Thus, the groundwater level should be adopted as the loading/unloading response for the calculation of LURR.

In this study, we attempt to devise an approach to determine the stress state of source media by applying the LURR method to identify the difference between groundwater levels during the loading and unloading periods induced by earth tides. The method for differentiation of loading and unloading phases follows the methodology of Coulomb failure stress change given by Yin et al. [29]. To show the validity of the approach, the groundwater level in the wells near the epicenters of the 2008 Ms = 8.0 Wenchuan, 2014 Ms = 6.5 Ludian, and 2019 Ms = 6.0 Changning earthquakes in the Sichuan and Yunnan provinces, southwest China are chosen as examples. The text of this paper is divided into seven sections. In addition to the introduction section, Section 2 and Section 3 show the method for LURR calculation and the preprocessing of the water level data; Section 4 and Section 5 display the application of the method to seismic data and some obtained results, and in Section 6 and Section 7 we would like to give some discussion and conclusions.

2. Methods

Our approach is founded on the premise that the increased volume in rock resulted from the cracks generated by the establishment of the criticality of an earthquake can change groundwater levels in the nearby wells. According to the studies of Scholz et al. [30], when tectonic stress reaches a relatively high level, the rock will become dilatant (Figure 1) due to the occurrence of cracks, i.e., volume magnifies relative to the elastic variations [28]. Unlike the linear elastic stage (Figure 2a), the increased volume in this stage can alter water level in the near field wells by changing permeability and pressure. On the other hand, in terms of the definition of the Kaiser Effect [31], the fracturing of rock during the loading and unloading periods is different under cyclical tide stress. Since the cracks tend to occur in the loading rather than the unloading [32], the difference between water levels during the loading and unloading periods is observed (Figure 2b).

Thus, the average values of groundwater levels recorded in corresponding periods are adopted to compute LURR. Specifically,

$$Y=\frac{\left({{\displaystyle \sum }}_{i=1}^{N+}{H}_i\right)/\left(N+\right)}{\left({{\displaystyle \sum }}_{i=1}^{N-}{H}_i\right)/\left(N-\right)}$$

(1)

where H_i is groundwater level at the i-th record, and N+ or N− represent respectively the numbers of data points during the loading and unloading periods. The time window used to calculate LURR values usually involves multiple load-unload cycles to reduce violent fluctuations in the curves.

Since cracks in rocks are usually caused by shear stress [33], in our model, the Coulomb failure stress (CFS) induced by earth tides in the tectonically preferred slip direction on certain slip surface is used to determine stress change of source rock, and the Coulomb failure criterion is adopted to differentiate the loading and unloading phases (as Yin et al. [34] did in their LURR calculation using Benioff strain of small earthquakes as the loading/unloading response).

In the investigation of earthquake preparation, two stresses should be taken into consideration: the tidal and tectonic stresses [35]. The effective shear stress (τ_e) in earth crust is written as the common effect of the tectonic-effective shear stress ( ${\stackrel{\rightharpoonup }{\tau }}_e^T$) and tidal-effective shear stress ( ${\stackrel{\rightharpoonup }{\tau }}_e^t$):

$${\tau }_e={\stackrel{\rightharpoonup }{\tau }}_e^T+{\stackrel{\rightharpoonup }{\tau }}_e^t.$$

(2)

This shear stress is also called the CFS which is usually expressed as:

$$CFS=\tau -f\sigma ,$$

(3)

where $f$, $\tau$, and $\sigma$ denote, respectively, the inner frictional coefficient, shear stress, and normal stress on the stress plane [36].

The orientation of tectonic-effective shear stress, which can be regarded as the principal stress orientation of CFS [37], is determined by the slip vector ( $\stackrel{\rightharpoonup }{u}$) on the fault surface of the ensuing mainshock [38]. ${\stackrel{\rightharpoonup }{\tau }}_e^T$ and ${\stackrel{\rightharpoonup }{\tau }}_e^t$ can then be defined as (Figure 3a):

$${\stackrel{\rightharpoonup }{\tau }}_e^T={\stackrel{\rightharpoonup }{\tau }}^T-f{\sigma }^T\frac{\stackrel{\rightharpoonup }{u}}{\left|\stackrel{\rightharpoonup }{u}\right|},$$

(4a)

$${\stackrel{\rightharpoonup }{\tau }}_e^t={\stackrel{\rightharpoonup }{\tau }}^t-f{\sigma }^t\frac{\stackrel{\rightharpoonup }{u}}{\left|\stackrel{\rightharpoonup }{u}\right|}.$$

(4b)

On the other hand, because the loading rate of tidal stress is much greater than the tectonic stress [35,39], i.e., $\frac{\Delta {\tau }_e^T}{\Delta t^T}«\frac{\Delta {\tau }_e^t}{{\Delta t}^t}$, the change of effective shear stress on the principal stress orientation can therefore be expressed as:

$$\Delta {\tau }_e={\stackrel{\rightharpoonup }{\mathit{\tau }}}_e^t·\frac{\stackrel{\rightharpoonup }{u}}{\left|\stackrel{\rightharpoonup }{u}\right|}.$$

(5)

Generally, the sign of $\Delta {\tau }_e$, which is decided by the angle of θ (Figure 3b), determines the loading or unloading stressed by earth tides. When θ < π/2, loading, and θ > π/2, unloading.

If the rate of CFS is taken into consideration, the loading and unloading can then be decided by the sign of g which is defined as:

$$g=\frac{d}{dt}({\stackrel{\rightharpoonup }{\mathit{\tau }}}_e^t·\frac{\stackrel{\rightharpoonup }{u}}{\left|\stackrel{\rightharpoonup }{u}\right|}).$$

(6)

The method for the calculation of tidal-induced stress is similar to the methodology proposed by Yin et al. [34]. The detailed algorithm has been listed by Yu et al. [40], in which the elastic deformation in crust induced by a celestial body is written as 6 differential equations of first-order [41], and the stress components on any section can be derived by using the Runge-Kutta method [42].

3. Data Preparation

The observation of groundwater level is subject to the monitoring instruments, artificial noise, and environmental influences, all of which cause errors of some kinds. Hence, the water level data are preprocessed by following procedures to enhance the quality of the LURR calculation.

• To remove data extrema from the water level sequences. For the rock media, the magnitude of increased volume result from the cracks, in terms of the studies of Brace et al. [25], is no more than 2.0 times the elastic volume variations. Thus, the data points whose values exceed twice the average amplitude of the water level sequence are removed.
• To linearly interpolate values to the missing data in the water level sequences.
• To perform 12.4~26 h (mainly including the M2 and O1 tidal data) Butterworth bandpass filtering on the data to remove signals unrelated to the tidal process.

4. Application to Seismic Data

As a retrospective study, we apply the LURR approach to the groundwater level in the wells near the epicenters of the Ms ≥ 6.0 earthquakes in the Sichuan-Yunnan region of southwest China. The study region is located geographically in the southeast corner of the Qinghai-Tibet plateau (Figure 4). The Qinghai–Tibet plateau is derived from the sustained collision and compression between the Indian and Eurasian plates. Because of the strong tectonic movement of the Qinghai–Tibet plateau, the Yunnan–Sichuan region is seismic active, in which more than 30 Ms ≥ 7.0 earthquakes have been recorded.

The China Earthquake Administration (CEA) has monitored groundwater in this region for many years for seismic hazard evaluation. Nevertheless, few earthquake cases can be used in this study because our method is limited by the spatial distribution of observation wells. The near-field water level observations with a relatively clear tidal process are indispensable. There are only three earthquakes with suitable water level observations that can be applied to calculate LURR: the 12 May 2008 Ms = 8.0 Wenchuan, 3 August 2014 Ms = 6.5 Ludian, and 17 June 2019 Ms = 6.0 Changning earthquakes.

The Wenchuan earthquake occurred in the middle and northeastern parts of the Longmenshan (LMS) fault zone which lies in the topographic boundary between the eastern margin of the Tibetan Plateau and Sichuan Basin. The LMS fault zone consists of the low-angle transpressional faults which extend more than 300 km from southwest to northeast. The seismogenic fault of the Ludian earthquake is the Baogunao-Xiaohe fault, with about 10 km rupture. The Baogunao-Xiaohe fault is a northwestern trending sinistral strike-slip fault which is located in the east part of the border between Sichuan and Yunnan provinces. The Changning earthquake occurred in the basin-mountain junction zone on the southern margin of the Sichuan Basin. The epicenter is located on the Changning-Shuanghe anticline, with the seismogenic structure mainly caused by the NE-SW compression.

The groundwater levels at six wells (Figure 4) were selected for the calculation of LURR. As regards the Ludian earthquake (LDEQ), the epicentral distances of the ZYLJ and ZT wells are 35.8 and 42.8 km. Larger epicentral distances of wells are observed for the Wenchuan earthquake (WCEQ) and the Changning earthquake (CNEQ). More specifically, regarding the WCEQ, the epicentral distance for the QL well is ~76 km, while for the CNEQ the epicentral distances of wells ranging from 71 to 106 km.

Table 1. Detailed information of the well water level adopted for LURR calculation.

Event	Well Name	Epicentral Distance (km)	Well Location (Lat/Lon)	Start Time	Depth (m)	Sampling Frequency
LDEQ	ZYLJ	35.8	27.38/103.58	20120524	114	1-h
	ZT	42.8	27.45/103.58	20100101	350	1-h
CNEQ	LZ	85	27.38/103.58	20130524	300	1-h
	RCHJ	106	27.45/103.58	20140101	251	1-h
	NX	71	28.98/104.93	20120601	105	1-h
WCEQ	QL	76	30.32/103.28	20070630	175	1-h

Note: The groundwater level observation at the NX well was stopped on 1 April 2019.

lists the detailed information of the observation wells. The raw data of groundwater levels were displayed in Figure 5. The wells were installed from 2007 to 2014, whose depth goes from 100 to 400 m. The resolution of the probes is shown in

Table 2. The instruments and their resolution of the geophysical observations in Figure 4.

No.	Name of Station	Observation Type	Instrument	Resolution
1	Zhaoyangleju	Groundwater level	TDL-15	0.001 m
2	Zhaotong	Groundwater level	LN-3A	0.001 m
3	Luzhou	Groundwater level	SWY-Ⅱ	0.001 m
4	Rongchanghuajiang	Groundwater level	SWY-Ⅱ	0.001 m
5	Nanxi	Groundwater level	LN-3A	0.001 m
6	Qionglai	Groundwater level	SWY-Ⅱ	0.001 m
7	Shuifu	Horizontal pendulum tiltmeter	SSQ-2	0.0005″
8	Zhaotong	borehole strain	RZB-2	5 × 10−9
9	Wenchuan	Horizontal pendulum tiltmeter	JB-1	0.002″
10	Guza	borehole strain	YRY-4	1 × 10−10

. Figure 6 shows some examples of the preprocessed water level data, which are prepared for LURR calculation. It is clear that the earth tide can be observed in the water level fluctuations.

In order to show the noise components in the water level due to non-tectonic factors, we perform tidal analysis by using the Fast Fourier Transform (FFT). The frequency of the different harmonics and their amplitudes are shown in Figure 7, in which the semidiurnal and diurnal tide components can be easily picked up. The figure is calculated by using the tidal harmonic analysis method suggested by Pawlowicz et al. [43].

Figure 8 shows the time series of LURR produced using the groundwater level as data input, with the time window of 60 days at a sliding step of 30 days. Detailed source models and slip distributions (used for computing CFS) for the three mainshocks are adopted from the global CMT solution (http://www.globalcmt.org/, accessed on 10 January 2021), and the internal friction coefficient is 0.4 [34].

5. Results and Interpretation

Comparing the LURR time series derived from the groundwater level at different wells, we find that significant anomalies can be observed prior to the ensuing large earthquakes. The LURR time series before the LDEQ began to increase in late 2012 and reached maximum values in early 2014 (Figure 8a), while the LURR peaked in early 2019 and in late 2007 for the CNEQ and WCEQ, respectively (Figure 8b,c). Note that the annual changes could be observed in the groundwater levels of LZ and RCHJ wells (Figure 5d,e) and their corresponding LURR time series (Figure 8b), which may be caused by the precipitation. To clearly show the variations in LURR due to the dilatancy, we calculated the influences of annual changes on the LURR (Figure 8b). The annual changes in water level can be obtained (removed) by using the Slipping Fourier Analysis method proposed by Zhao and Liu [44]. By comparing LURR derived from the water level in which the annual variation is removed with the LURR shown in Figure 8b, the errors of LURR caused by the annual variations could be determined (Figure 8b). Even if the errors are deducted from the LURR peaks, marked anomalies in the LURR time series before the earthquake are also clearly observed. The results could be supported by previous works of Yin et al. [34], which illustrated that the LURR evolved around 1.0 for many years until 1~2 years before the mainshocks when, the LURR reached maximum values (about 1.19, 1.08, and 1.1 for the three earthquakes, respectively), and subsequently return to a low level several months prior to the earthquakes.

Figure 8 also shows the monthly rainfall recharges in the areas where the observation stations are located. During the period of anomalous increase of LURR, the rainfall did not change (increase or decrease) significantly. Therefore, the pre-seismic changes in the LURR time series should not be caused by the effect of rainfall recharge on the loading/unloading of the aquifer.

On the other hand, there still are noticeable differences between the LURR curves: [1] for the event with multiple observation wells, the magnitude of LURR anomalies decreases with the distance from the epicenter: the shorter the distance, the higher the peak. The LURR peak values prior to the LDEQ are 1.19 and 1.16 for the ZYLJ and ZT wells respectively, while the corresponding distances from the epicenter are 35.8 and 42.8 km (Figure 8a). Similar scenarios can be observed in the LURR time series before the CNEQ (Figure 8b). [2] For some large events in the same neighborhoods, significant pre-seismic anomalies may also be detected not just to the target earthquake. As shown in Figure 8c, in addition to the WCEQ in 2008, the LURR time series still yielded significant anomalies using the water level in the QL well before the 2013 Ms = 7.0 Lushan earthquake (LSEQ). The WCEQ and LSEQ ruptured respectively the northeastern and southern segments of the Longmenshan fault zone (Figure 4). The distance between the two epicenters is just 86 km, and the QL well is located in the middle of them, with distances of about 76 km to the WCEQ and 28 km to the LSEQ. Due to the occurrence of the WCEQ, the increased stress on the southern part of the Longmenshan fault [45] induced a relatively high value of LURR until the LSEQ.

It is quite interesting that the stress state change in the crust identified by LURR may also be detected by the near field geodetic time series in the same time periods (Figure 9). We found that the borehole strain (BHS) of ZT and GZ, and the horizontal pendulum tiltmeter (HPT) of SF and WC are about 49, 81, 57, and 56 km from the epicenters of the LDEQ, LSEQ, CNEQ, and WCEQ (see Figure 4). The marked compression in the ZT and GZ BHS in late 2012 (Figure 9a,b), the obviously uplift changes in the SF HPT during 2018~2019 (Figure 9c), and the remarkable deflection in the WC HPT in 2007 (Figure 9d) provide potential evidence of volume variations for the pre-seismic dilatancy, and anomalous changes of LURR produced using near field well water level. Note that the earth tide is clearly recorded in these geodetic time series (Figure 10).

6. Discussion

The LURR is a usable and convenient parameter for exploring crustal stress state. The LURR presented in this paper is evaluated based on the water level changes derived from the increased volume of cracks within the source rock. Since the tide stress is far less than the tectonic stress, cracks can only be triggered, not created, by tidal stress. When tectonic stress is at a low level (e.g., the elastic phase in Figure 1), it is hard to induce any cracks by tidal stress. There were no noticeable differences between water levels during the loading and unloading periods, so that Y ≈ 1.0 (Figure 2a). When tectonic stress reaches a relatively high level, the source rock is very sensitive to any tiny stress change, and cracks may therefore be triggered by the CFS result from tidal stress. Combined with the studies of the Kaiser Effect, we know that cracks tend to be triggered by the increased CFS, i.e., under cyclical tide shear stress, more cracks occur during the loading than the unloading, resulting in the detectable difference in the groundwater level (Figure 2b), and Y > 1.0.

Our studies show that, for an earthquake with M ≥ 6.0, the LURR has very significant prediction efficiency when the epicentral distance is less than 100 km, and the efficiency decreases obviously when it exceeds this distance. The volume changes associated with the water level (LURR) fluctuation could be found in the geodetic time series in the same neighborhoods. We note that the ZT BHS began to show marked compression at the end of 2012 (Figure 9a) when the LURR values produced using the groundwater level of ZYLJ and ZT were greater than 1.02 (Figure 8a). Similarly, the EW and NS components of the SF HPT (Figure 9c) together with the LURR values derived from the LZ, NX, and RCHJ wells began to increase in 2018 (Figure 8b), while the time of LURR anomalies before the 2008 WCEQ (Figure 8c) could be approved by the uplift changes in the EW and NS components of the WC HPT (Figure 9d). A suitable explanation for the changes in the geodetic and LURR time series at almost the same time is the pre-seismic dilatancy [28]. The physical basis of the dilatancy model has been discussed by Wawersik and Brace [46]. Cracks increase the void volume and change the stress state of the source rock so that the anomalous changes in the geodetic time series were observed, and then with the influx of water, the increase of LURR values is observed. The co-seismic water level fluctuations were clearly observed in the wells selected for LURR evaluation (Figure 5). This could make it possible to consolidate the idea that these boreholes could be representative of stress state variation of the aquifer. However, these changes cannot be detected in the time series of LURR, because the co-seismic water level changes have been removed in data preprocessing.

In practice using the approach, when LURR is greater than a statistically critical value ($Y^{\mathrm{c}}$), does tidal triggering of cracks occur? By comparing the curves displayed in Figure 8 and Figure 9, we can see that when the geodetic time series is running in a stable mode, the LURR value is basically less than 1.02 (indicated by the horizontal dashed lines shown in Figure 8). When the geodetic time series anomalously changes, the LURR is usually greater than 1.02. The results correlate well with the evolution of the rock constitutive curve (Figure 1). During the establishment of dilatancy, the difference between the responses to loading and unloading is augmented gradually. That is, with increasing external load, the LURR value should gradually increase to the maximum. Thus, the LURR anomalies can last months to years and peaked 2~8 months before the earthquake when the load reaches the limitation of P_max (Figure 1). Then, due to the reduction of the response difference between the loading and unloading, the value of LURR is dropped prior to the mainshocks. Thus, 1.02 might be set as the optimally critical value of LURR for testing events in the Sichuan-Yunnan region.

Compared with the raw data of groundwater level (Figure 5), the pre-seismic changes identified by using the LURR technique look very significant in the entire sequence, which manifests clearer warning signals for testing final events. The LURR generally fluctuated around 1.0, while the maximum value is about 1.1 (Figure 8). The difference between them is 0.1. On the other hand, the critical LURR value is 1.02, with a difference of 0.02 from the standard value. That is, the variation of the maximum LURR value is five times that of the critical value. The change makes the LURR an effective precursory parameter of earthquakes.

Finally, the calculation of LURR requires a relatively precise assessment of the direction of the tectonic-effective shear stress which depends on the location and focal mechanism of the ensuing large event [47]. For a given earthquake fault, if it is driven toward failure, its surrounding areas should be loaded to an anomalously high level by the tectonic stress. The high stress state in areas surrounding the seismogenic fault is manifested by the sensitive changes of water level in the near field wells response to the cyclically tidal loading before the mainshocks. Because the water level data are specially differentiated into the loading and unloading phases by using the effective shear change induced by the cyclical tide stress in the tectonically preferred slip direction on certain slip surfaces, rock dilatancy or development of cracks can be effectively identified, and the LURR anomalies for testing events in the same neighborhoods are therefore presented. Thus, where enough knowledge of active tectonics and stress setting exists, the LURR approach, combined with the high-quality water level observations, can be tuned toward detecting the potential of future large earthquakes.

7. Conclusions

We find that by calculating the ratio between the groundwater levels during the loading and unloading periods induced by tidal stress, marked changes might be detected prior to the occurrence of a large earthquake. For an earthquake with Ms ≥ 6.0 in the Sichuan-Yunnan region, the critical LURR value is about 1.02 and the anomalies can last several months to years. The changes could be attributed to the pre-seismic dilatancy, with the volume changes in the earth crust recorded by the geodetic time series in the same neighborhoods. A combination of the pre-seismic changes in the LURR and geodetic time series would represent less ambiguous alarms for detecting the upcoming large earthquakes. At present time, the application of this method for real-time prediction of the earthquake must be approached with caution because pre-determined source models are needed. Nevertheless, the attempts shown in this paper make us systematically search for the time and location of various kinds of earthquakes by using the near-field water level, if sufficient knowledge of regional tectonic stress setting is known. More importantly, our approach elucidates an important feature of groundwater level that could be combined with the physically feasible models for seismic hazard evaluation.