The Longmenshan Fault, positioned critically between the Tibetan Plateau and the Sichuan Basin, extends impressively over a length exceeding 500 km and spans approximately 70 km in width, featuring a northeast-southwest strike (Figure 1). The fault zone is subject to continuous compression from the Tibetan Plateau in the northwest, resulting in highly active geological activity (Zhang, 2008). Throughout the Late Quaternary period, the fault zone’s activity has exhibited a gradual intensification from north to south. The fault zone is developed within a metamorphic heterogeneous rock body, characterized by high rupture intensity, thereby facilitating energy accumulation and predisposing the area to the occurrence of powerful earthquakes. Since the 1960s, “Y”-shaped fault zone (F1, F2, and F3 in Figure 1) has experienced a total of seven earthquakes with magnitudes of 7.0 or greater in Sichuan Province, establishing it as the most active region for strong earthquakes in western mainland China (Bai et al., 2019). The 2008 Wenchuan Ms8.0 earthquake has generated large ruptures of up to 300 km in length beneath the surface, occurring within an exceptionally brief time frame. The central rupture zone has been observed to span approximately 240 km (Zhang et al., 2008).

Emerging from the southernmost segment of the Longmenshan fault zone, the Anninghe fault constitutes an additional region prone to frequent seismic activity. This fault, extending in a north-south direction, is predominantly characterized by sinistral strike-slip faults. It spans approximately 170 km in length and exhibits a complex hierarchical structure (Yi et al., 2004). Occupying a pivotal role, the Anninghe fault has been the site of a series of earthquakes (He and Ikeda, 2007).

Following the statistical relationship between magnitude and distance as detailed in the “China Seismic Code,” the study selected earthquakes at a certain distance from the monitoring wells, near the Longmenshan-Anninghe Faults (98–107°N, 24–34°E). Specifically, earthquakes of magnitude Ms5.0-6.0 were chosen within a 250 km seismic distance from the wells, earthquakes of magnitude Ms6.0-7.0 within 300 km and those of Ms7.0 or higher within 500 km (Table 1). This provides reasonable assurance that all monitoring wells will be within the range of seismic static strain (Shi et al., 2013). The earthquake time frame selection criteria are based on recent earthquake occurrences (2017–2023) and the inclusion of a wide range of magnitudes, particularly focusing on earthquakes above magnitude 5 and 6 to ensure methodological compatibility. Moreover, the groundwater level time series is selected within the same 2017-2023 timeframe to correspond with the seismic events under investigation. The surface wave magnitude (Ms) was chosen as the magnitude type for this paper, which was determined by measuring Rayleigh wave amplitudes at periods of approximately 20 s. Unlike Local magnitude (ML), which is calibrated for regional distances (<600 km) and specific to local crustal structures. Moment magnitude (Mw), based on the seismic moment tensor solution, provides the most complete physical description of earthquake size by considering fault parameters. However, Ms remains the standard scale in China, particularly effective for shallow earthquakes (depth <70 km) and historical catalog comparisons. A total of eight groundwater monitoring wells, depicted in Figure 1, were chosen for this study to ensure a better correspondence with the earthquakes. Some wellbore and stratigraphical are shown in Figure 2. Monitoring wells were deployed by the China Earthquake Administration (CEA) and were equipped with LN-3 and ZKGD3000-N groundwater level detectors, recording at a frequency of one measurement per minute, 1 mm resolution and 0.2% F.S. The original groundwater level curves and the difference curves are illustrated in Figure 3. Difference curves are calculated using first-order differences, which helps in highlighting changes in the rate of change of ground water level. The key information regarding the monitoring wells and their related features are outlined in Table 2.

The Molchan diagram method offers a quantitative approach to assess the correspondence between groundwater level changes and subsequent events (Molchan, 1990; Zechar and Jordan, 2008). It involves establishing various differential thresholds to calculate the Abnormal time period occupancy rate τ and the Miss rate v can be calculated. These values are then plotted as τ-v step lines within the Molchan diagram, also known as the Molchan test line. The position of the step lines determines the strength of the correlation between groundwater level changes and subsequent earthquakes (Molchan, 1990; Zechar and Jordan, 2008). Molchan diagram requires the assessment of probability Gain and significance, and the equations involved and the significance of the parameters are as follows (Zechar and Jordan, 2008): G a i n = h τ = 1 − v τ (1) B h | N , τ = C h N τ h 1 − τ N − h (2) where h is the number of hits: the number of earthquakes that successfully landed in the alarm region; v is the miss rate: The ratio of earthquakes not falling within the alarm region to the total number of earthquakes; τ is the abnormal time period occupancy rate: the ratio of the anomalous time horizon of the groundwater level to the total; B is the cumulative binomial distribution, which is used to test for statistical significance; and N is the number of random hits. Gain is determined by the combination of v and τ, and the length of the time period does not affect the results. The closer the Molchan test line is to the line of greater probability gain, the better its overall prediction. For convenience, we define the normalized area to the right of the Molchan test line as the pre-response index (PRI), that is, the potential of groundwater level to reflect subsequent earthquakes. The PRI range is 0–1, and the closer it is to 1, the stronger the correlation between groundwater level changes and subsequent events. The PRI value changes with the position of the Molchan test line (Sun et al., 2017).

In this study, a criterion is established to differentiate between high and low PRI. This criterion is based on setting a threshold value that reflects the absolute magnitude of differential values. Differences that exceed the threshold value are identified as anomalous. The alarm area will be set within a certain period of time after anomalous. An earthquake is deemed to have been successfully hit if it occurs both within the alarm area. Conversely, a hit is regarded as unsuccessful if the earthquake happens outside of the alarm area. To facilitate understanding, let’s consider the example of well C-48 for a 94-day period from 1 July 2018, to 1 October 2018 (Figure 4). We will set the alarm area as 10 days.

a) When the threshold line is set to 0.5, only the value of July 3 is determined to be anomalous. At this time, according to Equations 1, 2, τ = (1 + 10)/94 = 0.12. The 7th earthquake does not fall within the alarm area, thereby v = 1/1 = 1.

In use, there is no need to manually select a threshold. The Molchan method automatically traverses the cycle from the differential water values maximum to the minimum value, we can obtain multiple sets of v and τ corresponding to the different thresholds. Molchan diagram takes into account the combined results of all thresholds and avoids subjectivity in identifying anomalies.

Wavelet coherence analysis quantifies the correlation between groundwater level and theoretical tidal series by measuring their temporal relationship (Grinsted et al., 2004; Song et al., 2023; Yang and McCoy, 2023). It identifies resonance periods through phase-shifted arrows in highlighted regions, revealing the degree of correlation between two time series X and Y: R 2 s = S s − 1 W X Y s 2 S s − 1 W x s 2 * S s − 1 W y s 2 (3) where W ^X and W ^Y are discrete wavelet transforms, W ^XY is the cross wavelet transform of X and Y, S is the smoothing window, and R ² is the coherence coefficient. R ² is ranging from 0 to 1, with values close to 1 indicating that groundwater levels and tides vary in a high correlated manner. The wavelet coherence coefficient resembles the correlation coefficient in the traditional sense, and it can be understood as a localized correlation coefficient within the frequency space. Simply enter two time series groundwater level and theoretical tide with the same resolution, and R ² between them will be calculated according to Equation 3. The code based on MATLAB is already publicly available for download (Grinsted et al., 2004).

Based on the tidal effect of groundwater level, the Leaky Aquifer Model is used to invert the specific leakage. First, tidal analysis of the water level data facilitates the determination of both observed and theoretical values for various tidal sub-waves’ parameters. In this analysis, two key parameters are the amplitude ratio, which is the observed amplitude divided by the theoretical amplitude, and the phase shift, which represents the difference between the observed phase and the theoretical phase. Both parameters are essential for understanding the tidal analysis. By constructing response models for different well-aquifer systems, the phase shift and amplitude can be utilized to invert the hydraulic parameters of the aquifer.

In situations where the aquifer’s water primarily flows horizontally towards the borehole, a negative phase shift is observed. The radial flow model (Cooper et al., 1965; Hsieh et al., 1987) can be employed to invert the permeability coefficient under these conditions. However, in more realistic scenarios where the aquifer interacts with surrounding rocks through hydraulic processes such as leakage, the phase shift tends to exhibit a leading behavior. The leaky aquifer model can be utilized to derive the specific leakage (σ), expressed as σ = k'/b', where k' and b' represent the permeability coefficient and thickness of the aquitard (Gu et al., 2024; Wang et al., 2018). Significantly, this model also accounts for scenarios where flow within the aquifer is purely radial (k' = 0). The specific leakage (σ) serves as an indicator of the aquitard’s vertical water transport capacity. The theoretical equations of leaky aquifer model are as follows: T ∂ 2 h ∂ r 2 + 1 r ∂ h ∂ r − K ′ b ′ h = S ∂ h ∂ t − B k u ρ g ∂ ε ∂ t (4) where T (m²/s) and S are the transmissivity and storage coefficient of the aquifer, respectively, r is the lateral distance from the well, k and k' are the permeability coefficients of the aquifer and the aquitard, respectively, b and b' are the thicknesses of the two, and B and K _u are the skempton’s coefficient and the undrained bulk modulus, respectively. Equation 4 has the analytical solution as: A = a b s i ω S i ω S + k ′ / b ′ ξ (5) η = arg   i ω S i ω S + k ′ / b ′ ξ (6) Where ξ = 1 + r c r ω 2 i ω r ω K 0 β r ω 2 T β K 1 β r ω (7) β = k ′ T b ′ + i ω S T (8) where A and η are the tidal parameter amplitude and phase shift, respectively, and r _c and r _ω are the case pipe radius and well filter pipe radius, respectively.

In use, we need to input the known parameters: A, η, r _c , r _ω and Bku, then the specific leakage can be computed by combining Equations 5–8 in MATLAB using open-source code (Zhang et al., 2024).

Random forest regression (RFR) was used for sensitivity analysis and was able to quantify the potential contribution of multiple factors to PRI (Borup et al., 2023; Rigatti, 2017). It uses the bootstrap resampling technique to generate a new set of training samples by repeated random sampling of n samples from the original training sample set T. Each independently sampled training sample is used to train a tree, and the n decision trees generated from the sample set are computed in parallel to select the optimal result, which improves the model’s generalization ability. The Gini index is used to complete the establishment of the regression tree, the smaller the Gini index, the better the decision tree division (Breiman, 2001). Assuming that the sample T contains k classes, the Gini coefficient can be expressed as: G i n i T = 1 − ∑ k = 1 n p k 2 (9) where p _k denotes the probability that the sample belongs to the kth class. The smaller the Gini index, the smaller the uncertainty will be and will be more useful for feature testing. After normalizing all the factor series and inputting them into the RFR at the same time as the PRI, the contribution of each type of factor to the PRI will be output. The code is based on MATLAB's own function Treebagger, and the resluts of Equation 9 will automatically calculate (https://ww2.mathworks.cn/help/stats/treebagger.html).

Some studies have reported changes in groundwater levels prior to certain earthquakes (Wang and Manga, 2021), but these observations are not universal or consistent. The mechanisms of such changes remain poorly understood, and currently there is no reliable way to use groundwater level variations as earthquake precursors. Therefore, it is necessary to statistically screen valid anomalies from a large number of suspected anomalies and establish the correspondence between groundwater level anomalies and subsequent earthquakes by the Molchan method. Anomalies are defined in this paper as differences exceeding a certain threshold. In the Molchan method, the threshold is adaptively selected (more details see Section 3.1). The premise of the method: a certain number of samples are needed, the samples include groundwater level and seismic data, and the longer the groundwater level series, the better. In the groundwater level time series, we can identify suspected anomalies, but this may not actually be the case (Figure 3). After accumulating a certain number of suspected anomalies, they are then matched with seismic data occurring in the vicinity. The more samples involved in the calculation, the more accurate the statistical results will be, ultimately revealing the correlation between groundwater level changes and subsequent earthquake events.

Figure 3 shows that pre-processing groundwater level data using differential values effectively highlights anomalies in groundwater levels. These anomalies are mainly marked by notably high values that deviate significantly from the surrounding data points, particularly in the period preceding the earthquake event (Differential values in green circles). Wells C-43, C-42, and C-47 exhibit more anomalies prior to the earthquake, potentially providing clearer indications of seismic activity. In contrast, wells C-05, C-03, and C-48 show no significant groundwater level changes, and the remaining two wells display only general fluctuations. However, subjectivity is not a discriminating criterion, and we will use the Molchan diagram to a further test.

Molchan method unifies the characteristics of groundwater level changes before earthquakes, which provides a quantitative indicator that we defined as pre-response index (PRI). PRI represents the area to the right of the Molchan test line, ranging from 0 to 1.The closer it is to 1, the stronger the correlation between groundwater level changes and subsequent events (Lai et al., 2021; Sun et al., 2017). Figure 5 shows the results of the PRI for eight wells under a 30-day alarm region period. The v-t test lines for wells C-41, C-47, and C-48 are closer to the Gain line of Gain = 2, and grey areas are higher and more significant, indicating a relatively higher PRI. The v-t test lines for wells C-05 are the lowest, falling well below 0.5, indicating that the PRI is relatively lower.

Furthermore, to enhance the precision and reliability of the Molchan test, an assessment of PRI was conducted across various alarm regions, limited to a maximum of 180 days. Figure 6 presents the PRI of eight monitoring wells within this timeframe. Distinct variations in PRI were observed among the wells. The average PRI over a 6-month period ranged from a minimum of 0.35 to a maximum of 0.71, indicating a considerable disparity. Wells C-41, C-43, and C-47 have a good performance in general (PRI of more than 2/3), while the PRI of wells C-03 and C-05 is relatively low (PRI close to 1/3), C-42, C-46 and C-48 show average performance (between 2/3 and 1/3).

Despite the similarity in tectonic units and the wells’ location in fractured media aquifers, as well as the comparable distances between the selected earthquakes and the wells, the groundwater level PRI exhibits noticeable variations. These variations provide a unique opportunity to identify factors that govern groundwater level changes before earthquakes.

The aquifer characteristics, particularly the confinement, play a crucial role in determining the wells’ responsiveness to seismic stress. Furthermore, hydraulic parameters significantly influence the magnitude of water level response, thereby affecting the characteristics of groundwater level changes. Tidal effect and seismic static stresses induce certain disturbances within the aquifer and can be considered as “typical representatives” of aquifer response to external stresses. Groundwater levels that exhibit favorable responses to both types of stresses are more likely to demonstrate satisfactory correlations with subsequent earthquakes. Additionally, the distribution of epicenters in relation to the stress propagation path may also influence PRI, serving as an important factor that warrants consideration. In this section, we selected the degree of tidal action experienced by the wells, the co-seismic response magnitude, the hydraulic properties of the aquifer, and the orientation of the epicenter as factors to correlate with PRI. This comprehensive analysis aims to shed light on the factors contributing to the difference of PRI among the monitored wells. The PRI serves as a reliable indicator of the correlation between groundwater level changes and subsequent earthquake events.

This analysis involved multiple steps. First, we evaluated the co-seismic and tidal response coefficients by basic statistics and wavelet coherence analysis, since the role of tides is potentially significant and needs to be treated in the frequency domain to highlight the correlation between the time series. In addition, spectral analysis and leaky aquifer model were utilized to assess the aquifer confinement and hydraulic parameter, both of which are suitable for dealing with periodic signals similar to tidal action. Based on these evaluations, a comparative analysis was performed to assess whether the conditions were responsible for the observed differences in PRI. Secondly, we quantitatively analyzed the impact of the epicenter’s distribution location on the PRI of the wells’ water levels. This analysis aimed to identify any correlations between the spatial distribution of seismic events and the PRI. Lastly, we employed the RFR method to integrate the analyzed factors with the sensitivity analysis of the factors controlling the differences in PRI. This is because it is suitable for multiple series to be analyzed simultaneously with strong robust-ness. This combined analysis aimed to identify the factors that are more likely to contribute to the variations in PRI. Ultimately, our goal was to identify the key factors that significantly impact the variations in PRI.

Hydraulic parameters and seismic distribution can significantly affect the PRI of water level, and the degree of water level response to stress can also reflect the PRI, while the confinement of the aquifer in wells may not have a significant effect on the PRI. In order to identify the relative influence of various factors on the PRI of water level, the sensitivity analysis of these four types of factors was conducted using the RFR. This method effectively integrates both quantitative and qualitative results, and is easy to be combined with the previous analysis. Here, weights in confinement to unconfinement are given as 1, 0.5 and 0, and other relevant factors need to be normalized as well. Then by inputting the series corresponding to each type of factor and PRI, RFR can automatically generate results.

Table 5 presents the weight proportions obtained from the RFR. Among the selected factors, distribution of the earthquakes stands out as the predominant control, exerting a substantial influence on the discrepancies observed in PRI. The hydraulic parameters also contribute to the variations in water level PRI, but with a weightage that is only slightly lower than the distribution of the earthquakes. In contrast, the degree of water level response to external stress (Tidal effects and Co-seismic responses) and the confinement conditions of the aquifer exhibits comparatively weaker impacts on PRI, as indicated by their lower weightages, indicating their relatively lesser significance in determining the PRI of water level.

In our selected cluster of wells, C-46, C-42, and C-48 exhibiting moderate PRI in comparison to the remaining wells. These three wells are situated at the intersection of the Longmenshan-Anninghe faults, thus are significantly affected by the distribution of earthquakes because seismic signals are more likely to propagate near the rupture zone, but our screening conditions for seismic events are not oriented. If earthquakes do not occur near the target faults, potential earthquake precursor information is just as easily lost.

Furthermore, differences between wells C-41, C-43, and C-47, characterized by high PRI, and wells C-05 and C-03, identified as having poor PRI, can largely be attributed to the differences in hydraulic parameters of the aquifer, since most of which are located far from faults intersections. The specific leakages differ by 2–3 orders of magnitude, as confirmed in the previous section.

The sensitivity of water level to external stresses (Tidal effects and Co-seismic responses) also offers some insights into PRI, although it is as a secondary criterion. In this regard, the contrasting performance of excellently responding well C-43 compared to wells C-05 and C-03 underscores the potential influence of stress sensitivity on PRI, and can be a useful reflection of the overall expected performance of the well. In our selected wells, the aquifer confinement on PRI shows some randomness, and it may show correlation on a larger scale, so our results conclude that confinement is not enough of a key cue for high or low PRI.

Furthermore, additional factors such as well depth, borehole radius, and the height of the water column within the well should be considered, as they may potentially contribute to variations in PRI. However, their impact on water level PRI is relatively minor. For instance, the slight differences in PRI between wells C-41 and C-43, which have having similar PRI, or among wells C-46, C-48, and C-42, where the difference in PRI varies only slightly within a range of 0.05, illustrate the limited impact of these factors. Although these factors can be relevant in certain scenarios, their overall contribution to altering water level PRI remains secondary in comparison to the dominant factors previously discussed. Hence, the factors affecting the disparities in PRI, ranked in descending order of sensitivity, encompass 1) distribution of the earthquakes in relation to well locations, 2) hydraulic parameters and 3) the sensitivity of water level to external stress responses.