The Pollino 2012 seismic sequence: clues from continuous radon monitoring

The 2012 Pollino (Calabria, Italy) seismic sequence, culminating in the Mw 5.2 earthquake of 25 October 2012, is investigated, exploiting data collected during a long-term continuous radon monitoring experiment performed in the epicentral area from late 2011 to the end of 2014. We analyse data collected both using a phenomenological approach based on quantitative evidence and a purely numerical analysis including the following: (i) correlation and cross-correlation investigations; (ii) an original approach aimed at limiting the impact of meteorological parameters variations on the interpretation of measured radon levels; (iii) a change point analysis; (iv) the implementation of an original detection algorithm aimed at highlighting the connections between radon emission variations and major seismic events occurrence. Results from both approaches suggest that radon monitoring stations can be subject to massive site effects, especially regarding rainfall, making data interpretation harder. The availability of long-term continuous measurements is crucial to precisely assess those effects. Nevertheless, statistical analysis shows a viable approach for quantitatively relating radon emanation variations to seismic energy release. Although much work is still needed to make radon time series analysis a robust complement to traditional seismological tools, this work has identified a characteristic variation in radon exhalation during the preparation process of large earthquakes.

station the reduction effect of rainfall on radon observations seems less marked, but it is still present. Figures 3a, 3b, 3c and 3d show selected fall-winter and spring-summer periods for MMNG, respectively. In this case, though the reduction of radon emission with precipitation is still present (Fig. 3a), heavy rain events cannot be clearly separated from radon concentration peaks, being sometimes overlaid (yellow rectangles) (Figs. 3b,c and d).
3 Analysis of radon timeseries 5 In the following we try both to assess the impact of meteorological parameters on radon signals on a quantitative basis and to outline an original approach aimed to remove (or at least mitigate) the effects of meteorological events on the detected timeseries. Our goal is to maximize the informative power of radon emanation variations potentially related to a variation in seismic energy release.
Even though the effects of meteorological conditions on temporal radon timeseries have been investigated for the last fifty 10 years by means of different approaches and methodologies (Singh et al., 1988;Zmazek et al., 2003;Piersanti et al., 2015), a clear assessment and a solid interpretation has not been univocally established yet.
For the following analyses, we decided to use only radon timeseries from station MMN, since it was the only one installed before the main events of the sequence (M W 4.3 on May 2012 and M W 5.2 on October 2012), corresponding to the major changes in cumulative seismic moment release rate (Fig. 4c). From data collected in the time window from April 2012 to 15 December 2012 (Fig. 5a), that includes the two major seismic events, we note that in correspondence of these two changepoints the radon emanation increased a few days before the seismic events. Both the average amplitude and duration of such increases appear to scale with the magnitude of the corresponding earthquakes, as highlighted in the two yellow rectangles of Fig. 5a. The apparent discontinuity in the radon increase just after the M W 5.2 seismic event is likely to be associated with a major precipitation episode right after the earthquake occurrence. observed values continue to increase up to about 1600 Bq m −3 for more than 30 days after the earthquake (except, as described above, for the first days of November, when a major precipitation event flattened down radon levels).

Correlation and cross-correlation analysis
In order to quantitatively assess the phenomenological evidences described above by means of numerically objective procedures, we perform a series of statistical evaluations on our dataset. Figure 4 shows the whole timeseries employed in the statistical analysis filtered with a 14-days moving average. All the moving averages employed in our computations are evaluated backwards (i.e., average at day d i , employs only the previous (d i -14) days). Figure 4a represents the radon concentration 30 (black line) and rainfall (red line), Fig. 4b shows temperature (black line) and pressure (red line). Figure 4c shows the cumulative seismic moment release (black line) with the seismic moment release (red line).
Since the Pearson coefficient reflects mainly a linear relationship between variables, we estimated the correlation between variables using both the Pearson coefficient (Hollander et al., 2014) and a non-parametric correlation coefficient (Kendall, 1970). The two approaches yield virtually identical results, so we show here only the classical Pearson analysis. We performed both a correlation analysis between radon and environmental parameters and a cross-correlation analysis between radon, meteorological parameters and seismicity. All analyses look for a linear relationship between two variables but the cross-correlation 5 considers it as a function of the time-offset of one relative to the other. Formally cross-correlation function reads (i.e., Chatfield,

2004):
where N is the series length, u t and y t are the two time series, u and y are their sample means, and k is the lag. Differently from Pearson linear correlation, the cross-correlation coefficient is not normalized a-priori: in order to grant compatibility with 10 the previous analyses, we normalized the cross-correlation coefficient here so that it varies between -1 and 1 and set the lag range between -40 and 40 days.
We decided to exclude rainfall from this analysis since, differently from other meteorological variables, it is intrinsically characterized by a strongly discontinuous, spike-like behavior being the majority of the sampling times characterized by a null value. In fact, during the time window of our most relevant analyses, we have null rain values ranging from 65 % to 75 % of 15 the sampling intervals (to compare for instance with less than 10 % of days with null seismic moment release). This makes correlation and cross-correlation analysis inadequate approaches to evaluate the relationship between radon concentration and rainfall.
The results regarding the correlation analysis in terms of Pearson coefficient are summarized in Table 1. For all considered cases, we report both global cumulative value (G) corresponding to the entire acquisition window (2012-2014) and separate 20 results for each year (2012-2013-2014). The Pearson coefficient ρ shows a significant level of negative correlation only between radon concentration and temperature with values ranging from −0.6 to −0.2. The value of the Pearson coefficient for pressure, even though coherent both in sign and in magnitude for each time window, is nevertheless statistically compatible with zero.
Within the cross-correlation analysis, whose results are shown in Fig. 6, we include also the seismic moment release M 0 , since for this physical variable a lagged approach is able to consider also a causal relationship in addition to an instantaneous 25 feedback among variables (Box and Jenkins, 1976;Piersanti et al., 2015). Figure 6 is arranged in nine panels: from left to right the cross-correlation between radon and temperature, pressure and seismic moment release, respectively, are presented, while the rows represent 2012, 2013 and 2014 time-windows. No sharp and isolated peak is observed in Fig. 6, indicating that no clear cross-correlation scenario can be deduced from this analysis. Nevertheless, we can confirm the correlation pattern the Pearson coefficient ρ for which the probability of obtaining a cross-correlation greater than or equal to ρ for uncorrelated data is equal to 1 (Chatfield, 2004) and is represented in Fig. 6 by the grey lines. The cross-correlation function between radon concentration and seismic activity shows a significative positive peak during 2012 (when the major seismic events occurred), with a maximum value of 0.5 in correspondence of a 21 days delay forward of radon concentration (Fig. 6, panel upper right).
Of course the relationship between variations of radon emanation and seismotectonic processes would be better assessed if 5 we would be able to remove, or at least reduce, the bias of meteorological parameters on the radon measured concentration.
To this aim, we implement an empirical correction procedure for temperature, pressure and precipitation variations. Basically, given an observed radon concentration value Rn obs , taken at time t when a temperature T , an atmospheric pressure P and a precipitation level R have been registered, we define a corresponding meteorological-corrected concentration Rn cor as: where C P , C R and C T are positive correction factors obtained as a simple linear interpolation from the minimum detected 10 values of T, P and R in a selected time window where {C P = C R = C T = 1} (that is to say there is no-correction), to the maximum detected values in the selected time window where {C P = C Pmax ; C R = C Rmax ; C T = C Tmax }. The optimal value of C Pmax , C Rmax and C Tmax can be obtained by maximizing the cross-correlation function for the selected time window (of course a time window including a significant seismic activity must be selected). We want to note that the subscript max above stands for maximum magnitude of the correction, not for maximum absolute value of the correction parameter C i . Indeed, if 15 the correction factor corresponding to the maximum value of a given meteorological parameter C i is > 1, it means a negative correlation between radon and that parameter, the opposite if the correction factor C i is < 1. Since it is reasonable to consider the possible connection between radon concentration variations and seismotectonic processes as dependent from the seismic source-observer distance (Dobrovolsky et al., 1979), we have implemented in the correction procedure also the possibility of weighting for the epicentral distance (Hauksson and Goddard, 1981;Einarsson et al., 2008). Again, given an earthquake with 20 seismic moment M 0 obs occurred to an epicentral distance r from station MMN, we consider a corresponding distance-weighted value M 0wgt : where w is a positive weighting factor (w=0 means no correction for epicentral distance).
In Fig. 7 we show the effects of our correction procedure on the cross-correlation function. The extrapolation of the optimal values for the correction parameters C Pmax , C Rmax , C Tmax and w was performed by means of the MINUIT package (James, from April 2012 to January 2013. As can be seen from Fig. 7, the proposed correction procedure significantly increases crosscorrelation peaks for both time windows (indicated as tw-1 corrected and tw-2 corrected). Notably, the increase is greater for the larger time window where a lower (but still significant) peak cross-correlation value was obtained, while the time lag of the peak remains completely unchanged after the correction, indicating that the variation of radon intensity seems to follow the variation in seismic moment release. In Table 2 the correction coefficient values maximizing the cross-correlation peak in the 5 two time windows tw-1 and tw-2 are reported. From the tabulated values we note that: i) the correction values for the rainfall lie in both cases at the top of the searching domain (C Rmax =10 for tw-1 and C Rmax =9 for tw-2), i.e., rainfall is strongly anticorrelated with radon emanation, confirming the phenomenological analysis in previous Sect. 2; ii) the correction values for the temperature are always greater than 1, confirming that for MMN station temperature is anti-correlated with radon emanation (see above in this same section); iii) the correction values for the pressure oscillate about C Pmax =1, confirming the lack of a 10 clear correlation regime between pressure and radon emanation for this station.

Change point analysis and detection algorithm
The problem of detecting changes in timeseries is well known in climate literature: the definition and identification of discontinuous steps, or change points, may be subjective and it also depends on the form of the trend one expects between changes.
Several methods have been implemented to solve the change point problem both for short and long climatic timeseries. We refer 15 the readers to Reeves et al. (2007), in which the literature about the change points methods is widely reviewed and discussed.
We applied to the measured radon intensity timeseries an algorithm developed in the realm of Earth's climate system studies in order to calculate, by means of a bayesian approach, the posterior probability of multiple change points in a generic climatic timeseries (Bayesian Change Point algorithm, (Ruggieri, 2013), BCP hereinafter). Once the algorithm has identified an arbitrary number of change points in our timeseries, whose maximum is an input parameter of the algorithm (k max = 6 20 in the following) , our primary interest is to verify if the detected change points in the radon timeseries are consistent with corresponding changes in cumulative seismic moment release rate (i.e., major earthquakes).
Applying the BCP algorithm to the whole MMN timeseries, we obtain an indication of most likely two change points that are potentially associable with the two largest events of the sequence. Figure  be useful to get further insight into the relationship between radon and seismicity, employing the same BCP algorithm on the the cumulative seismic moment release timeseries, in order to check the possibility of finding significant variations in seismic moment release other from the trivial ones (i.e. coincident with a major seismic event). The result is shown in Fig. 9: in this case the rates changes, that are clearly visible a priori, are all found by the algorithm with a probability near to 1. While the 2 nd and the 4 th change points clearly identify the two earthquakes, the 1 st and the 3 rd change points seem, instead, to identify 5 the beginning of a preparatory phase of the two events. The first occurs on February 20, 2012 (1 st red spike in Fig. 9) and the third on August 18, 2012 (3 st red spike in Fig. 9). We note that the temporal difference (about 70 days) between each of these two change points and the change points estimated by the BCP algorithm for the MMN timeseries (the two blue dashed vertical lines in Fig. 9) is comparable. In this respect, radon concentration variations could be sensitive to the internal processes taking place during the preparatory phase of an earthquake.

10
We point out the fact that a standard change point analysis uses always the whole timeseries, since to identify a change point at a time t i the algorithm processes also data at t > t i . This is a limitation because the algorithm cannot be employed for predictive purposes. To overcome these limitations and most of all to extend the range of our investigations, we implemented an original detection algorithm that potentially could be used in real time analyses. A schematic flow chart of the algorithm is shown in Fig. 10.   to that observed by means of change point analysis.
Therefore, both the cross-correlation analysis and the change point analysis, as well as the application of our detection algorithm, indicate that a physical relation between the variation of soil radon emanation and seismic moment release exists.While change point and detection algorithm both succeed in finding some useful radon signal before the variation in seismic moment release, the cross-correlation investigations seem to behold the radon signature after the seismic moment release variation. Relying on the change point analysis and detection algorithm, we have verified if also the cross-correlation analysis is compatible with a radon signal preceding the seismic moment release signal. To investigate this possibility, we have repeated the procedure described in Subsect. 3.1. In this case we limit the search domain to positive lag values (i.e. radon signal preceding moment release signal), in order to verify if a suitable solution can be found also in this case. As Fig. 12 highlights, such a solution exists and, comparing Figs. 7 and 12, it is evident that it is only marginally less significant with respect to the best one.

5
Remarkably, as a confirmation of the previous findings, the correction coefficients associated with this solution (see Table 3) are consistent with these found in Subsect. 3.1. They indicate for radon observations at MMN station a strong anti-correlation with respect to precipitation (C Rmax =9.3 for tw-1 and C Rmax =10.0 for tw-2), a clear anti-correlation with temperature and the lack of a clear correlation with respect to pressure variations. 10 We have performed a detailed analysis of the temporal variations of radon emanations from late 2011 to 2014 in a seismically active area during a seismic sequence that culminated at the end of 2012 with a M W 5.2 event. We exploited several different approaches to carry out our investigations. Namely: i) phenomenological analysis; ii) correlation and cross-correlation investigations; iii) empirical correction of the meteorological parameters effect on radon timeseries and its impact on crosscorrelation; iv) change point analysis; v) detection algorithm. 15 We can split the main results of our work in two classes: a) those concerning the impact of meteorological parameters variation on the observed radon timeseries and b) those concerning the existence of a physical connection between the observed radon timeseries and the seismic moment release temporal variations. Converging indications coming from both classes represent an important outcome of our work. Regarding class a), we have indications that, in the investigated setting, soil radon emanation is strongly anti-correlated with precipitation and weakly anti-correlated with temperature, while we do not get sig-20 nificant and univocal evidence of correlation (positive or negative) with pressure variations. In this context, approaches i), ii)

Conclusive remarks
and iii) give remarkably consistent indications and we see as particularly significant the agreement between the strength of the correlation evidenced by i) and ii) and the magnitude of the corresponding correction factor found with iii). These results, when compared with previous findings, confirm that the environmental impact on radon observations is strongly site dependent. The correlation between radon variations and temperature is, in this sense, a clear example: many works found it positive, 25 as several others (including ours), negative. This observation suggests that a specific characterization is needed for each station, when implementing an observational network (see, for example, the dependence on the varying soil characteristics as porosity, permeability, and pre-rain moisture state). Regarding class b), all our analyses univocally indicate the existence of a non-accidental correlation between the temporal evolution of soil radon emanation and seismic moment release. The primary output of approach ii) suggests that the radon signal follows the seismic moment variation, while approaches i), iv) and v) 30 indicate that it is possible to retrieve the radon signal also before the seismic moment variation. Remarkably, we have found that even if approach ii) gives as primary result a shifted forward temporal correlation, nevertheless, also the solution with the radon signal preceding the seismic moment variation is acceptable at a barely lower significance level.   Table 2.  Table 2. Correction coefficients for temperature (CT max ), pressure (CP max ), rainfall (CR max ) and epicentral distance (w) maximizing the cross-correlation function (CC) in time-windows tw-1 (Jan2012-Jan2013) and tw-2 (Apr2012-Jan2013).