Implication of Radon Monitoring for Earthquake Surveillance Using Statistical Techniques: A Case Study of Wenchuan Earthquake

The seismotectonically induced changes in groundwater radon (Rn) are considered to be strong imputes for the surveillance of imminent major earthquakes. Groundwater facilitates the migration of soil gases as a result of tectonic stresses. In this regard, a radon time series is statistically analysed to identify the radon anomalies possibly induced by Wenchuan earthquake. The statistical analysis mainly involves the deterministic analysis of the Rn data and residual Rn analysis using a criterion ð x ± 2σÞ of anomaly selection having a confidence interval of 95%. The deterministic analysis reveals that the Rn time series follows a persistent trend ð0:5 ≤H ≤ 1Þ which confirms the absence of a chaotic regime. On the other hand, the residual Rn shows a notable upsurge straddling the time of the Wenchuan earthquake in the form of preand post earthquake changes at monitoring stations having RE/RD ≤ 0:3. The residual Rn level passes the anomaly selection criterion ð x ± 2σÞ and is declared as a tectonically induced Rn anomaly. Contrary to this, the response of distant monitoring stations (RE/RD > 0:3) to this particular earthquake further validates the link between Rn and earthquake activity. In a nutshell, the present study highlights the potential implications of earthquake-induced radon anomalies for earthquake prediction research.

Among them, the hydrogeochemical response of earthquake activity has been extensively studied for the last few decades [24]. The hydrogeochemical studies mainly include the monitoring in water level changes, Rn variability in either air, soil, or groundwater, an emanation of soil gases (CO 2 , CH 4 , and He), and other constituents such as Cl − , SO 2− 4 , and HCO − 3 are also used as a forerunner to access the earthquake-induced changes [7,10,11,[24][25][26][27]. However, the use of Rn is globally preferred in comparison with other precursors due to its relatively long half-life and ease of predictability [14].
In the past few years, numerous physical models were proposed to explain earthquake precursors [1-3, 20, 28, 29]. For example, Pulinets and Ouzounov [3] proposed a Lithospheric-Atmospheric-Ionospheric-Coupling (LAIC) model that explains the synergy between earthquake activity and its precursors. According to the LAIC model, the relative movement of tectonic blocks leads to the generation of tectonic stresses. Furthermore, Koike et al. [30] performed a laboratory experiment to verify the growth of Rn emissions as a result of compressional stresses. These tectonic stresses are responsible for opening and closing of microfractures/cracks that facilitate the migration of Rn from the crust to the upper surface by means of carrier fluids such as groundwater or other soil gases like CO 2 , CH 4 , or N 2 in a wide range of geological settings [20]. In particular, Cicerone et al. [2] reported a correlation between the concentration of Rn and the regional variation of stress/strain. In addition, numerous studies highlight the presence of notable Rn anomalies in groundwater around the time of earthquake activity [19,[31][32][33].
In connection to the above, the Tashkent earthquake (M 5.3; 1966) provided the first evidence of earthquake-induced Rn anomalies in groundwater and gave hope to research on earthquake precursors [34]. Several years' analysis of the Rn data shows a multifold increase in the Rn levels around the time of this particular earthquake, and the same patterns also repeated for its aftershocks over a short time scale [35]. Such an analogous behaviour of Rn, earthquake, and its aftershocks gives evidence of earthquake-induced Rn anomalies. Later, a successful prediction was made for the Haicheng earthquake (M 7.3; 1975) based on multiple datasets along with Rn as discussed by [36]. Similarly, Hirotaka et al. [37] observed a gradual increase in the Rn count two weeks prior to the Nagano Prefecture earthquake (M 6.8; 1984) 65 km away from its epicentre. Afterwards, Igarashi and Wakita [5] analysed the Rn anomalies associated with earthquake events occurring from January 1984 to December 1988 in Japan and its surrounding areas. Prior to correlating Rn with earthquakes, the Bayesian statistic was successfully applied to remove background variations along with a criterion for anomaly selection (2σ), event magnitude (M ≥ 6), and epicentral distance (D ≤ 100) from the monitoring site. Twenty Rn anomalies satisfied the anomaly selection criteria and exhibited notable pre-and post earthquake changes, respectively. The Kobe earthquake (M 7.2; 1995) was also preceded by a marked groundwater Rn anomaly at a 17 m deep observation well 30 km away from the epicentre [38]. Presently, numerous studies have reported the groundwater Rn anomalies and suggest a more detailed analysis of the Rn anomalies for a more reliable prediction of major devastating earthquakes [11,17,20,27,33,35].
In this regard, the present study provides a detailed temporal analysis of the Rn fluctuations in groundwater straddling the time of the major (Mw 7.9) and shallow (depth ≈ 19 km) Wenchuan earthquake. The Wenchuan earthquake occurred on May 12, 2008, along the Longmenshan fault (31.0°N, 103.4°E) in Sichuan Province, China [39,40]. It is the most destructive earthquake in China since the 1976 Tangshan earthquake [41]. This devastating event was responsible for thousands of casualties and massive destruction on a regional scale. The Wenchuan earthquake also triggered a large number of geohazards such as landslides, slope collapses, and debris flow that seriously damaged counties in Sichuan Province, China [42,43]. Given the devastation caused by this event, it is imperative to increase knowledge on earthquake precursors in this region. Therefore, we have analysed the groundwater Rn variations at near and distant locations. It is important to mention that numerous researches also analysed the fluctuations of different ground and space-based parameters; however, the systematic analysis of radon is absent. Therefore, the present study will help to establish a synergy between numerous parameters in the context of earthquake forecasting.

Seismotectonic Settings of the Investigated Region
Seismically, China is ranked among highly active regions of the world with extreme interplate seismicity. The famous plate tectonic theory is unable to explain this intraplate seismicity of China, which state the occurrence of earthquakes within a narrow plate boundary zone [44]. The historical catalogue of China reveals the occurrence of more than 1000 (M > 6) earthquakes since 23 A.D., with thirteen being classified as a major catastrophic (M > 8) event. A classic example of a devastating historical earthquake in China is the Huaxian earthquake (1556) which killed approximately 830,000 people and is considered the deadliest earthquake in human history [45]. The Tangshan earthquake (M 7.8; 1976) is the best-known example of a modern age event in China in terms of loss of life (killed~250,000 people and injured millions) and property [46]. In particular, the seismicity of western China is mainly influenced by the Indo-Eurasian collision which is responsible for the generation of the E-W trending fault system [47,48].
The key tectonic feature in this region is the Longmenshan fault (LMSF) that runs parallel (around 500 km long and 40-50 km wide) to eastern Tibet and Sichuan basin in the northeast and southwest directions [12]. The LMSF zone is mainly dominated by the Wenchuan-Maowen fault, the Yingxiu-Beichuan fault, and the Guanxian-Jiangyou fault ( Figure 1) with several thrust faults resulting from the compression between the Yangtze craton and the eastern Tibetan Plateau [40]. Historically, only two notable earthquakes of moderate magnitude (M 6.5, 6.2) were recorded in 1657 and 1956 along the LMSF zone, and considered as seismically inactive, until 2008 [40,46], while Shi et al. [12] reported only one earthquake (Mw 6.1; 1989) on this fault. In 2008, the LMSF triggered a major (Mw 7.9) and shallow-focus (depth ≈ 19 km) Wenchuan earthquake (May 12, 2008) which gained the attention of the community in the capacity of reactivating the LMSF zone. The Wenchuan earthquake ruptured a 290 km long segment of the LMSF in total. This rupture is propagated 270 km unilaterally in the NE direction and 20 km in the SW direction with 80 km coseismic surface rupture. The Wenchuan earthquake resulted from the change of the LMSF dip angle (30°-50°SW to 60°-70°NE) and fault motion (SW thrusting motion to NE strike-slip motion) as provided commonly in various proposed slip models and summarized by [40].
The Wenchuan earthquake led to an estimated 69,225 causalities, 374,640 injuries, 17,939 missing, and approximately 5 million people left homeless. It is considered the second most destructive earthquake of this century after the great Sumatra earthquake in 2004 [41]. The Wenchuan 2 Geofluids earthquake was followed by more than 10,000 aftershocks of M ≥ 2:0 that lasted until Dec 31, 2008 [50]. This devastating earthquake was also responsible for the triggering of more than 15,000 geohazards (landslides, rockfall, and debris flow) which resulted in about 20,000 casualties [39]. The destruction caused by the event highlights the importance of earthquake precursory studies in seismically active regions to minimize the loss of life and property.

Methodology
The methodology adopted in the current study is divided mainly into two stages: (a) instrumental setup and (b) theoretical setup. The instrumental setup explains the strategy adopted for the data collection and a brief description of the instrument used for radon monitoring, while the theoretical setup explains about the statistical analysis of radon data for identification of radon anomalies possibly induced by earthquake activity. has started an organized, persistent, and systematic effort to testify the postulate of preearthquake anomalies in the context of earthquake forecasting by the China Earthquake Administration (CEA). This research programme was primarily focused on earthquake prediction and its implications to minimize the irrevocable destruction posed by major earthquakes. The CEA has mainly classified this research programme into the following classifications: imminent (weeks to days, even hours), short (months to weeks), medium (years), and long-term (decades) earthquake forecasting. In this regard, the CEA has installed seismological, geodetic/geoelectric/geomagnetic, ground fluid, and hydroseismological networks for monitoring potential earthquake precursors [42].
In particular, the CEA is operating a very dense radon monitoring network for the surveillance of earthquake activity. For earthquake forecasting research, Rn is more preferred due to its comparatively longer half-life and easy detectability [14]. In the current study, we have analysed the variations of continuous groundwater Rn concentration associated with the devastating Wenchuan earthquake ( L o n g m e n s h a n r e g i o n Q in g c h u a n fa u lt J i a n g y o u -g u a n g y u a n f a u l t J i a n g y o u -g u a n x i a n 3 Geofluids of the event epicentre. The details of the Rn monitoring stations regarding its location and distance from the event epicentre are presented in Table 1. The Rn data used in the current research was acquired by means of two instruments. For automatic continuous sampling, the mixture of escaped gas and water coming from well was passed through a degasser and gas-collecting device and then collected into a ZnS (Ag) detector system for Rn concentration measurement. This is a digitized method of measurement, and the observation equipment is mainly SD-3A with sampling interval of one hour and a measurement precision of 0.1 Bq/L. For the acquisition of the daily variation in Rn concentration, water from the well was sampled and degassed by bubbling degassing and then transported into an ion chamber or ZnS (Ag) detector, where the Rn concentration was measured by an ionization or scintillation method using FD-125 instrument. The measurement precision is also 0.1 Bq/L with a sampling interval of one reading per day. The further detail of instrumental setups used for Rn monitoring is provided by [51]. Precipitation is not monitored by the CEA for all the stations used in the current research, and therefore, such data is not available. To overcome this issue, the precipitation dataset of the Tropical Rainfall Measuring Mission (TRMM) satellite for the year 2008 was utilized in this paper. The TRMM dataset was considered as the best solution among the available sources of precipitation data for China. The negative correlation observed between groundwater Rn anomalies and rainfall for the whole year and the selected time frame is presented in Table 2.

Theoretical Setup
3.2.1. Fractal Analysis of Radon Data. The dynamics of the Rn time series show a very complex nonlinear temporal pattern normally characterized by nonstationary and multiscale features [52]. This chaotic regime of the time series is realized through diurnal, seasonal, multiyear, and decadal Rn cycles along with key influencing parameters [14,53,54]. Therefore, the Rn time series is subjected to fractal estimates to determine the degree of chaotic behaviour of Rn and intrinsic long-memory correlations, if any [55]. Besides this, the estimation of fractal elements for the Rn time series leads to further exploration of the underlying dynamics of physical systems such as seismic activity [56]. In this regard, the fractal quantity known as the Hurst exponent (H) is calculated for the Rn time series using the rescale-range (R/S) analysis. The estimation of H is based on the following relationships: where R is the range, S is the standard deviation, H is the Hurst exponent, and s is the number of entries in a group. where Here, the set of observations (N ∈ fN 1 , d with individual length = n/d. Afterwards, the Rn time series is categorized as antipersistent ð0 ≤ H ≤ 0:5Þ, random ðH ≈ 0:5Þ, and persistent ð0:5 ≤ H ≤ 1Þ based on the obtained value of H. In particular, the antipersistency means that low present values will probably be followed by high figure values and vice versa. Persistency exhibits that a long-lasting autocorrelation exists within the time series, which implies that high present values will probably be followed by high future values and vice versa. And random walk means that they are uncorrelated or do not possess long memory trend [14,56].

Identification of Radon
Anomalies. The inspection of fractal dynamics of Rn time series allows for the identification of anomalous Rn variations if they exist. For reliable identification of an earthquake anomaly, the monitoring station must lie within the Earthquake Preparation Zone (EPZ). The EPZ is defined as an area within which the premonitoring fluctuations associated with the tectonically induced impending earthquake can be observed. In [57], an empirical relationship is proposed for EPZ (R D ; km) based on the magnitude (M) of the earthquake event, given as where the epicentral distance (R E ; km) between the monitoring site and the event epicentre is calculated as where Here, α i , β i are the coordinates of the earthquake event, α j , β j are the coordinates of the Rn monitoring station, and R ≈ 6370 km is the earth's radius [58]. In an ideal scenario, only those events having R E /R D ≤ 1 are considered for earthquake forecasting studies. In our case, all the stations are lying within ðR D Þ as proposed by [57] and adopted in numerous studies [14,18,22].
In addition to all the above, the acquired groundwater Rn data is investigated for identification of anomalous periods possibly linked with this particular event. The recognized anomalous periods mainly lie around the time of the earthquake within a time window of ±3 months (Mar to Jul 4 Geofluids 2008). The rationale behind the selection of this flexible time window is to analyse the effect of whole aftershock sequence of the Wenchuan earthquake. The identified anomalous periods were further analysed using residual signal processing techniques to eliminate the regular filtering effect as indicated in numerous studies [14,22,59]. The residual Rn (dAðtÞ) is calculated via a relationship given as where AðtÞ is the daily average Rn concentration and RAðtÞ is the 7-day rolling average Rn concentration. Furthermore, we have also applied a statistical criterion ( x ± 2σ) of the anomaly selection on residual Rn having a confidence interval of 95% and is consistent with other studies [14,18,23,60,61].

Result and Discussion
The temporal variability of groundwater Rn was recorded from January to December 2008 near the LMSF zone at nine different locations in China. The data analysis reveals anomalous fluctuations of Rn at a few stations under normal meteorological conditions, which highlights the aspect of a tectonically induced Rn anomaly, while the inspection of the earthquake catalogue of this selected period owns the aspect of tectonically induced radon anomaly due to the Wenchuan earthquake that occurred on May 12, 2008, in Sichuan, China. Details of the earthquake are presented in Table 3.
In this regard, we performed a detailed statistical analysis of Rn concentrations, in order to verify the possible correlation of the Rn anomalies with this particular event. The Rn monitoring stations included in the current study are presented in Figure 2.
It includes 9 stations lying within the EPZ of the Wenchuan earthquake installed by the CEA within the proximity of the LMSF zone ( Table 1). The annual variation of Rn recorded at selected monitoring stations is presented in Figure 3.  (Figure 3(a)). This rise in the Rn level continues onward throughout the whole year. Likewise, other monitoring stations such as MXS, KDS, GS, and SPS also reveal an anomalous change in real Rn requiring a detailed study (Figures 3(b)-3(e)). An analogous change of the Rn level was also observed at the relatively distant monitoring stations with R E /R D > 0:5 (Table 1; Figures 3(h) and 3(i)). Based on the results of the preliminary investigation, the anomalous periods of the Rn concentration are subjected to a detailed analysis. This detailed analysis includes the advanced residual signal processing techniques which remove the regular filtering effects from data if any. Besides this, a statistical criterion of x ± 2σ is also applied to the residual Rn to further authenticate our results. The result of this detailed analysis is presented in two stages: (a) station is located very near to the event epicentre ðR E /R D ≤ 0:3Þ and (b) station is located far away from the event epicentre ðR E /R D > 0:3Þ for comparison purposes.
Prior to analysing the Rn data for earthquake forecasting research, we have determined the dynamics of the Rn time series using fractal dimensions such a Hurst exponent. The Hurst exponent reveals that the time series follows a persistent trend ð0:5 ≤ H ≤ 1Þ for all the Rn monitoring stations with insignificant fluctuations as presented in Figure 4. A positive autocorrelation is found to exist in the recorded Rn data. This suggests that the past trend of data is continued in the future and there is no existence of an irregular trend.

Geofluids
At the first stage, detail analysis of Rn is performed for the station located near the event epicentre as presented in Figures 5 and 6. It includes 7 monitoring stations located in the proximity of the LMSF zone having ðR E /R D ≤ 0:3Þ. Figure 5 illustrates the change of the residual Rn levels along with daily precipitation for MSS, MXS, KDS, and SPS stations. The temporal variation of the residual Rn concentration was observed from Mar 1 to June 30, 2018, along with daily precipitation records. The Rn levels seem to be within normal limits ð x ± 2σÞ until the occurrence of the Wenchuan earthquake on May 12, 2008. On May 12, 2008, the Rn levels breach the threshold of anomaly selection ð x ± 2σÞ and show a sudden upsurge from 0 to 2.4 Bq/L of the residual Rn and 13-18 Bq/L of daily Rn level (Figures 3(a)-5(a)). This particular change in the Rn level followed by the Wenchuan earthquake highlights the aspect of the post earthquake Rn anomaly, while the daily Rn variations recorded at the MXS stations show a gradual increase in Rn from April to June 2008 with insignificant precipitation (Figures 3(b)-5(b)).
A detailed inspection of Rn shows that the residual Rn level passed the anomaly selection criteria prior to the Wenchuan earthquake. This preearthquake Rn anomaly seems to be absent after event occurrence ( Figure 5(b)).
The KDS monitoring station shows unambiguous changes in the residual Rn levels around the time of this particular event (Figure 5(c)). Initially, this anomalous trend is observed 5 days prior to the Wenchuan earthquake on May 07, 2008, and continued until July 2008 showing post earthquake anomalies. It is important to mention that these post earthquake changes were possibly associated with the aftershock sequence of the Wenchuan earthquake [50]. Similarly, the SPS monitoring station having R E /R D ≈ 0:073 and R E ≈ 182 km shows a multifold increase in the daily Rn levels around the time of the Wenchuan earthquake (Figure 3(d)).   Table 1). The black dotted line is the time of the Wenchuan earthquake.

Geofluids
The average Rn value recorded at the SPS monitoring station ranges between 0.2 and 0.3 Bq/L, while the peak value observed one day prior to the Wenchuan earthquake is nearly 0.52 Bq/L. Additionally, the residual Rn level of the SPS monitoring station overpasses the anomaly selection threshold which further authenticates its linkage with this particular event.
In continuation, the temporal variability of the residual Rn concentration for GS, ZJS, and PZHS monitoring stations is presented in Figure 6. The GS monitoring station reveals an unambiguous increase in the Rn levels prior to event occurrence within a very short interval of time (4-6 days). The Rn level range is within a limit of 11-13 Bq/L from Jan 2008 to the end of April 2008. Afterwards, a sharp increasing trend of Rn is observed from the beginning of May 2008 and continued onward (Figure 3(e)). This anomalous change is subject to detailed analysis via residual signal processing technique for reliable identification of tectonically induced Rn anomaly. The results of the residual Rn show development of an Rn anomaly overpassing the anomaly selection threshold preceding the earthquake event ( Figure 6(a)). Similarly, the other monitoring stations presented anomalous patterns of the daily and the residual Rn levels highlighting the aspect of tectonically induced Rn anomalies under favourable conditions (Figures 3(f), 3(g), 6(b), and 6(c)). Inclusively, all the monitoring stations lying around the proximity of the LMSF zone show Rn variations with varying amplitude in connection with this particular event.
At the second stage, the temporal variability of Rn is analysed for distant Rn monitoring station with R E /R D > 0:3 (Figures 3-6). It includes the YQS and SYS monitoring stations having R E value 1542.9 km and 1563.2 km, respectively ( Table 1). The temporal variability of raw Rn levels recorded at these monitoring stations shows an increasing trend around the time of the Wenchuan earthquake, despite its distant location (Figures 3(h) and 3(i)). The YQS monitoring station shows multiple peaks in Rn levels with varying amplitudes prior to the event's occurrence. The highest peak (44.2 Bq/L) of Rn levels is observed on May 1-2, 2008, around 8 Geofluids 10-12 days prior to the Wenchuan earthquake, while the Rn levels remain within the normal limit (around mean value 39 Bq/L) for the rest of the period (Figure 3(h)). Likewise, the SYS station also shows an analogous change in the Rn levels around the time of this particular event (Figure 3(i)). These anomalous periods are further subjected to a detailed investigation as presented in Figures 6(a) and 6(b). The residual Rn value of the YQS station shows a significant increase on April 30, 2008, and overpasses the statistical criterion of the anomaly selection. This anomalous change in residual Rn is followed by the Wenchuan earthquake that occurred 11 days later, while during the rest of the period Rn levels were found to be within normal limits (Figures 3(h) and 7(a)). On the contrary, the temporal analysis of the residual Rn at the SYS station also depicts a notable rise in Rn levels almost 12-15 days prior to this devastating earthquake (Figure 7(b)).
Moreover, the comparative analysis of rainfall for the whole year and selected days further confirms the connection of Rn abnormality with this particular event ( Table 2).
It is important to mention here that the current results show an analogy with the earlier investigations of the preseismic process in connection with the Wenchuan earthquake [8,12,51,62]. For example, Shi et al. [12] studied the hydrological response of the Wenchuan earthquake at near and distant monitoring stations. They reported significant evidence of hydrological changes such as water level in response to the Wenchuan earthquake. Similarly, Ren et al. [51] also analysed the postseismic changes in Rn levels at different geochemical observation points using step variation curves. Their findings confirm the change in aquifer parameters as a response of dynamic loading, which may help in designing an optimum strategy for earthquake forecasting. In addition to this, Ye et al. [  9 Geofluids in soil and water along with water level and temperature fluctuations along the Longmenshan fault zone. They observed that both the Rn concentration and water level show positive and negative relationships in response to the Wenchuan earthquake. Besides this, Liu et al. [8] analysed the seismoionospheric parameter (GPS total electron content (TEC)) in connection with the May 12, 2008, Wenchuan earthquake. It is observed that the GPS-TEC shows an anomalous decrease above the forthcoming epicentral region and gives a possible validation of the LAIC model. In general, the results of the present study are found to be in accordance with earlier investigations and reveal that Rn is a promising tool for earthquake forecasting research.

Conclusions
The present study encompasses the analysis of the groundwater Rn response associated with the devastating Wenchuan earthquake (Mw 7.9) that occurred on May 12, 2008, in mainland China. On this subject, the daily variation of the Rn levels is recorded within the EPZ at multiple near and distant locations within the proximity of the LMSF zone that is responsible for the triggering of this particular event. Primarily, it is observed that the daily Rn levels show a notable change around the time of earthquake occurrence. The period of notable Rn variations is further subjected to the detailed analysis of Rn using statistical techniques for reliable identification of tectonically induced Rn anomalies.
The key findings of the current investigation are summarized below:    11 Geofluids (b) The deterministic analysis of the Rn time series suggests a positive autocorrelation with the absence of white noise in the dataset. The H exponent exhibits that the recorded data has a long memory correlation between past and future data points (c) The residual signal processing techniques were applied on observed periods of Rn enhancement for reliable identification of the Rn anomalies along with a statistical criterion of anomaly selection ð x ± 2σÞ. The response of the residual Rn variations associated with this tectonic activity appears in the form of possible pre-and post earthquake Rn anomalies with variable amplitudes. The findings of the current study are consistent with earlier investigations The analysis of the Rn time series shows a notable response possibly induced by tectonic origin that may effectively support the forecasting of impending catastrophic seismic events. Our results suggest that, on a global perspective, radon should be used as a potential seismic indicator for seismically active regions.

Data Availability
The data (groundwater radon) used to support the findings of this study belongs to China Earthquake Data Center and will be provided on demand from the relevant organization.

Conflicts of Interest
The authors declare no conflict of interest.