Stations CTTR and NRCA are equipped with a high sensitivity active radon monitoring instrument based on a Lucas cell15, with high efficiency (0.06 and 0.07 count min−1/Bq m−3, respectively) and minimum detectable concentration of 6 Bq m−3. The acquisition time is set to 115 min, and the stations acquire simultaneously local temperature values by means of a specific sensor. The daily values of external temperature, pressure, and precipitation are provided by the weather website https://www.ilmeteo.it as short term (12–24 h) forecasts. Figure S2 of the Supplementary material displays the available meteorological observables (temperature, pressure, rainfall) along with the radon time series collected at the two stations.
While the instrument at CTTR is located in a basement hosting the municipal archive of the city occasionally accessible only to technical staff (indoor installation type), the one at NRCA is co-located in a small shelter (shelter installation type) with a seismic station belonging to the INSN monitoring network. Systematic analyses of time series recorded by IRON stations indicated that indoor installations are noisier than shelter ones and show larger variability and stronger seasonal signal, likely connected with temperature fluctuations4.
This is confirmed by Fig. 2, showing for the two stations considered the distribution of radon concentration data through a density plot and a boxplot. The probability density function (Fig. 2a) estimated using Gaussian kernels16 indicates that the dataset recorded by the indoor station CTTR (blue) shows larger dispersion with respect to the NRCA one (red). The boxplots in the inset summarise the main parameters describing the statistics of the two datasets: CTTR data are more scattered, with a larger central value (median) and a smaller degree of right skewness compared to NRCA data. Figure 2b,c, where radon data from the two stations are grouped at month intervals, evidence that soil radon emission at station NRCA is maximum in the cold season (October-January) while the opposite holds true for station CTTR, where radon variations are more correlated with air temperature and more scattered during the whole year and markedly in the warm months. This seasonal periodicity has been also observed in laboratory tests17 and in several long term radon monitoring studies18,19,20.
Despite station CTTR has been recording almost continuously since July 2010, we selected here only the data collected since August 2012, when a new radon Lucas cell based detector was installed, for a total length of the time series of more than 7 years (2,367 days, 26,909 observations). Station NRCA was installed the day right after the occurrence of the 2016 Mw 6.0 Amatrice mainshock, so its data time series is less than 3-year long (1,042 days, 10,739 observations).
The time series of raw radon concentration data (Bq m−3) of the stations NRCA and CTTR are shown in Fig. 3a in red and blue, respectively. The longer time series of station CTTR shows a marked seasonal signal, less evident on NRCA one. Together with the different installation type, this is also due to the limited length of NRCA time series and by the presence of gaps in the measurements. A zoom on the time corresponding to the major seismic activity of the central Italy sequence (Fig. 3b) shows an apparently decreasing trend on CTTR data, possibly an effect of the seasonal behaviour. Among the various peaks in NRCA radon concentration, the two major ones are shared with CTTR time series and are therefore presumed to be significant. They precede the occurrence of EQ3, and are more pronounced (both in absolute value and relative to the background signal) in NRCA time series. The first radon concentration peak occurs 2 weeks before EQ3 and the second one occurs on the same day of EQ2.
Despite the wide literature on the subject18,19,20,21,22,23,24,25, a complete assessment of the relation between radon variations of tectonic origin and those due to atmospheric variables has not been established yet. In fact, the several studies addressing this question obtained results not univocal and often remarkably in contrast to each other. It is likely that the disagreement does not only reflect unreliability or weakness in the single analyses, but rather the strong site-dependency of the environmental effects, drastically influenced by the local features of the station site and by the installation typology.
A confirmation of this comes from the recent work by Siino et al. 2019, who carried on one the very few investigations in literature comparing as much as a dozen of concurrent radon time series widely distributed on a regional scale. The main outcome of this analysis, based on multivariate techniques, is that local site effects play a major role in modulating radon signal and, consequently, atmospheric correction approaches may not grant robust results.
Nevertheless, we make an attempt to employ an original empirical procedure14 aimed to limit the impact of environmental conditions on the interpretation of measured radon levels. The idea is to negatively weight the radon signal when it is associated to a peak of a meteorological parameter positively correlated with it, and to positively weight the radon signal if the meteorological parameter is negatively correlated with it.
To date, rainfall is the only meteorological parameter whose anti-correlation with radon data (due to the ‘ground sealing’ effect) is significant and widely observed. Conversely, the influence of temperature is strictly site-dependent, as it happens for the two stations considered here, which show opposite correlation with temperature variations (Figs. 2, 3). This is why at this stage of the analysis we chose to only take into account the effect of precipitation, leaving the treatment of the temperature- and pressure-induce bias to a future and more comprehensive study.
In Fig. 4, the raw and rainfall-corrected radon data, normalized according to their maximum value, are indicated in blue solid and dotted line, respectively. Green vertical lines indicate the time of occurrence of the four mainshocks and light-blue bars represent the daily amount of rainfall (mm). When the empirical correction is applied, the radon peak observed in raw time series on 15 October 2016 reduces of 50% in NRCA and becomes barely visible in CTTR, while the peak of 26 October still holds after the correction, meaning that it is presumably not associated with precipitation and rather to a deeper source common to the two stations.
Detection algorithm: approach
The visual inspection of the time series of radon concentration may provide some first order observations on radon emanation features, but the detection of more subtle changes in the time series for predictive purposes requires more advanced tools. We employ an original detection algorithm14 (DA hereinafter), which basically checks whether the radon daily average increases over a given threshold and if this condition persists over time.
The DA is aimed to highlight potential connections between radon emission variations and the occurrence of major seismic events. It might in principle be used in real time analyses since it explores the available data time series to forecast future seismic events issuing alarms on the basis of significant variations of radon emanation.
Piersanti et al. 2016 present in their Figure 10a flowchart of the algorithm operation: when the radon daily average exceeds the moving average (computed on p1 days) by a factor (p2) on the day i and the radon moving average successively increases by a factor (p3) for a time window (w = p4*p1 days), an alarm is issued at day i + w. If an earthquake of magnitude larger than 4 occurs within 40 days from the date of the alarm, the parameter p2 is increased by a factor p5 for a time window proportional to the earthquake moment. After that time, in case of no further M > 4 earthquake, the parameter p1 is restored to its original value.
Figure 5 illustrates how the tuning of the main input parameters (p1 and p2) affects the output of the algorithm. Red/blue circles indicate whether or not an alarm before EQ3 is issued running the DA on the radon time series of the stations NRCA (a) and CTTR (b). The number of combinations of the parameters that obtains an alarm are quite reduced in the CTTR case with respect to the NRCA one, meaning that the earthquake is more difficult to detect on that time series. This is likely due to the larger station-epicentre distance and to the lower SNR. When that distance is reduced, as for station CTTR relative to EQ1, the alarm before this earthquake is issued also for higher value of p2 (Fig. 5c).
In view of a real time application of the algorithm, another important factor we consider to select the set of input parameters for the DA is the number of false alarms delivered. Figure 5d,e, corresponding to NRCA and CTTR radon data, demonstrate that while the false alarms markedly depend on p2 (pardaily), they are only slightly affected by the value of p1 (durmedmob). Figure 5f reveals the higher number of false alarms issued for the CTTR time series compared to NRCA one, under the same value of p2. This confirms that clean and less noisy data are more suitable for the detection of tectonic signals.
The combination of parameters p1-p2 selected is highlighted in yellow and represents the tradeoff between the highest detection threshold p2 for which the DA finds an alarm before the earthquakes EQ1 and EQ3, and the minimum number of false alarms delivered. The whole set of input parameters is displayed in Supplementary Table S1.
Detection algorithm: results
Figure 6 shows the results of the DA applied to the time series of daily radon concentration (gray line) of NRCA (a) and CTTR stations (b). Green bars indicate the magnitude of the largest earthquake of the day occurred within a radius of 40 km from the station considered. The triangles represent all the issued alarms: in red the ones followed within 40 days by a Mw > 4 earthquake (yellow stars), in blue the remaining ones, deemed as ‘false’.
The DA applied to the NRCA time series successfully forecasts EQ3 and EQ4, along with a moderate-magnitude (Mw 4.4) earthquake preceding this last one. The first alarm is issued about 2 weeks before EQ3, a time advance consistent with the one characterizing the radon peak observed in NRCA and CTTR raw data (Fig. 3). We cannot rule out the hypothesis that the alarm points instead towards EQ2, occurred 4 days before EQ3.
As for the CTTR times series, despite the higher noise (Fig. 2), the DA succeeds in forecasting both EQ1, whose epicentre is closer than the one of EQ3 to this station, and EQ3, more far away but with larger magnitude. It also delivers 2 false alarms within the end of 2018, which rise to 6 when the DA analysis is extended to the whole time series of CTTR station (Fig. 7), for a total of less than one alarm/year.
These false alarms show a likely seasonal origin, being issued mainly in May/June or November/December, i.e. at the beginning or at the end of the ‘warm’ season, when radon emanation is observed to be larger at station CTTR. We performed a series of tests aimed to remove the annual periodicity of the radon data by subtracting the sinusoid which better fits the data, by differencing, by decomposition of the series (as a sum of trend, seasonality, and residual), by filtering out low frequencies (1 year) by means of the Fast Fourier Transform. After running the DA on the time series so processed, no alarm is issued at all (neither true or false). We interpret this result as an indication that the variability of radon is similar to the one of temperature, and since the signal is comparable to the noise, the detrending techniques end with cancelling out both.
Figure 8 is analogous to Fig. 6, with the difference that radon data are corrected for the effect of rain. The DA applied to the NRCA time series successfully forecasts the same events detected on the corresponding raw series, but issuing as many as 10 false alarms. The same behaviour is observed on the rainfall-corrected CTTR times series: the DA not only succeeds in detecting EQ1 and EQ3, but also delivers a further alarm for EQ4. The time advance of the alarm for EQ3 is about 2 weeks, consistent with that observed in Figs. 3 and 6. On the flip side, the number of false alarms throughout 3 years rises from 2 to 5, while on the whole 7-year-long time series its increase is of 6 times. We interpret this result as a consequence of the correction employed, which amplifies the values of radon concentration, leading the algorithm to detect anomalous radon variations more easily and more often. To get reliable results an ad-hoc tuning of the input parameters for the series corrected is unavoidable.
Change point analysis and cross-correlation
We applied to NRCA time series an algorithm originally developed in the realm of Earth’s climate system studies that calculates, by means of a Bayesian approach, the posterior probability of multiple change points (CP hereinafter) in a generic climatic time series27. The algorithm is able to identify an arbitrary number of CPs in time series, being the maximum number of allowed CPs and the minimum distance between adjacent CPs the input parameters of the algorithm (kmax and dmin, respectively).
Figure 9 shows the results of the Bayesian CP algorithm applied to NRCA time series with kmax = 3 and dmin = 2 with a number of sampled solutions fixed to 10,000; repeating the analysis with kmax and dmin in a reasonable range of values obtains similar conclusions. The algorithm finds a CP on 12 October 2016, i.e. preceding EQ3 of about 2 weeks, in agreement with what observed on raw data and what resulted from the DA analysis. The next 2 CPs found have higher a posteriori probability and likely correspond to the abrupt variations of radon emanation which reflect seasonal periodicity. This is observed also for the few CPs detected conducting the Bayesian analysis on the CTTR time series, not shown here for the sake of brevity.
To further check the delay time between variations in radon concentration and seismic activity, we cross-correlate the time series of radon data with the amount of seismic moment (M0) released every day by the earthquakes occurred on that day within a radius of 40 km from the station considered. M0 is obtained converting the total moment magnitude (Mw)28.
The time lagged correlation between the 14-day moving average of radon concentration and seismic moment release (Fig. 10) results maximum for a lag of -12 days. The negative offset suggests that radon is leading the interaction: correlation is maximized when M0 is pulled backward by 12 frames. The delay time is again in agreement with the time advance resulting from the DA and from the previous analyses.
The plot of the cross-correlation radon-seismic moment obtained for station CTTR (Supplementary Fig. S3) is maximized for a time lag of -8 days, which would be in agreement with the advance time of the alarm preceding the Mw6.5 mainshock (EQ3). But since the DA found two alarms (one before EQ1 and one before EQ3), we cannot attribute this delay to one or the other earthquake without the risk of over interpreting the result.