Introduction

Greek seismicity is among the highest worldwide due to the positions of the tectonic plates crashing at the east Mediterranean, giving options for radon studies in soil and groundwaters as an earthquake indicator (Fig. 1a). These experiments have been performed in Greece since the 80s and continue using the advanced monitoring technology of the time. The results contribute to the comprehensive worldwide database obtained over a half-century derived from experiments. Active faults have been investigated in Japan, Himalayan-Altai and Caucasus regions, Italy, Spain, and France [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. In addition, the seismic fault of St. Andreas in California gave the first results after long-term radon monitoring and the models between earthquake magnitude and the time passed from the radon signal, followed by models based on Himalayan studies [15,16,17,18,19,20,21,22,23,24,25].

Fig. 1
figure 1

a Greek seismicity, the marked area corresponds to the radon monitoring area b Location of radon stations (Gerakarou, Stivos and Scholari) and faulting studied [26]

The Greek peninsula has mostly been concerned with earthquake prediction since 2000, and integrated and continuous radon monitoring has detected abnormal changes in subsurface radon concentrations before earthquakes. The radon studies began after the catastrophic earthquake M = 6.5 occurred on June 20, 1978, in the Thessaloniki area and M = 6.6 on February 24, 1981, in the Aigion area. The Stivos fault at the Lagadas basin, near the Thessaloniki metropolitan area, has been periodically monitored from 1983 to 2002 [26, 27]. The Alkyonides and Aigion Faults at the Gulf of Corinth near the urban area of Athens and Patras were also occasionally studied from 1995 until today [28,29,30]. Another study in the Pyrgos area of West Greece near the Gulf of Corinth [31], on Lesvos Island and Crete [32], and an experiment in the Ioannina area in Northwest Greece [33] were also performed.

This work reports on data of radon variation in soil gas at the Stivos fault in the Lagadas basin associated with meteorological parameters and with data of seismic events to find radon signals predicting earthquakes. A time series analysis was performed at the continuous radon level data, temperature, and pressure. Radon signals after heavy rains were recorded, while radon anomalies to earthquakes were also detected. The results complete the events recorded at this fault since 1983, giving conclusions for modelling M with prediction time.

Instrumentation

The study was performed in the Thessaloniki area, North Greece, especially at the Lagadas Basin, monitoring the Stivos faulting in the Serbomacedonian zone in North Greece (Fig. 1b). Three stations for radon monitoring were placed at Gerakarou, Stivos and Scholari creating a triangle. Soil gas radon has been monitored since 1983–1985 and 1996–2000 using solid-state nuclear track detectors in specially made plastic devices 30 cm in height (see Fig. 2a). The device was set inside an outer plastic tube 1 m long with a 7 cm diameter [27]. From 1999 to 2002, silicon diode detectors in the BARASOL device 50 cm in height (ALGADE, France, see Fig. 2b) replaced the SSNTD, having a functional area of 450 mm2 of Si with 100 μm depth, a radon detection limit of 50 Bq m−3 and a saturation activity of 3 MBq m−3. The space 50 cm above the device was filled with Styrofoam material in granule form to isolate the devices thermally. Atmospheric parameters, temperature (°C), barometric pressure (mbar) and rain precipitation (mm) were recorded simultaneously by the sensors also included in probes [4, 7]. The devices were calibrated using a radon chamber of 110 l with a leakage rate 10–6 h−1, two orders of magnitude lower than the 7.5 10–4 h−1 radon decay constant, and a Ra-Rn source of 20 kBq [34].

Fig. 2
figure 2

Schematic diagram of a integrated radon monitoring device (using LR-115) buried in soil and b the continuous radon monitoring Barasol by ALGADE [26, 27]

Results and discussion

The results of radon concentrations in soil gas carried out with BARASOL devices for August 1999 to December 2002 at the Gerakarou, Stivos, and Sholari stations are shown in Fig. 3. The gap in the data collected from all radon stations (December 9, 1999–February 14, 2000) was due to software updates due to the new millennium. After 2000, technical problems occurred in the device installed at the Gerakarou station, causing the gaps observed. In the area of the M = 6.5 event in 1978, the radon levels at Gerakarou and Sholari were 25 ± 5 kBq.m−3 and at the Stivos 50 ± 15 kBq m−3 higher enough than the uncertainties involved in the measurements.

Fig. 3
figure 3

Radon level continuously measured using a BARASOL device associated with temperature, pressure, rainfall, and earthquake events at a Gerakarou, b Stivos and c Scholari station. The bold lines correspond to the daily moving average, while the oversized points to the rainfall and seismic events are discussed in the text

From August to December 1999, integrated measurements using SSNTD and continuous by Baracol device were in the field, presenting a good coherence between them (Fig. 4). The impact of rainfall and earthquakes is also clear, especially in the Stivos station, while no significant peaks and drops appeared at Scholari and Gerakarou (Fig. 3). A clear radon peak appeared at Stivos station due to rainfall in November 1999 and reached its maximum of 80 kBq.m−3 during two days after rains. The radon increases in the device after heavy rainfall due to increased barometric pressure inside the porous ground around the device that causes radon to flow through the device, which is not affected by the rain. After the radon peak, following the most interesting result, a radon signal reached its minimum of < 10 kBq.m−3 on November 26. Then it returned to normal levels on December 5, 1999—an inverse behaviour followed by the pressure giving a significant increase when radon was decreased dramatically. One week after the signal's end, an M = 4.5 earthquake occurred on December 12, 1999, predicting 36 days from the first sign and 7 days after the end of the radon anomaly.

Fig. 4
figure 4

Radon level continuously measured using a BARASOL device associated with temperature, pressure, rainfall, and earthquake events at Stivos during 1999

The radon peak due to rainfall was 2σ higher, while the radon drop detected before a seismic event was -3σ lower than the radon level standard deviation of 15 kBq.m−3. The effect of rainfall is shown in all stations, with radon peaks appearing during a few days after heavy rains. A rainfall radon peak + 2σ higher than the radon deviation of 5 kBq.m−3 appeared in Scholari station in April 2001 due to continuous weekly heavy rainfall. A peak was also detected in Stivos, but this time, the radon level was inside the radon deviation of 15 kBq.m−3. The above observation referred to any radon anomaly detected at Stivos station after 2000 instead of the clear anomalies detected during 1999. In the case of the Gerakarou station, no significant radon anomalies were also detected during 1999–2002 that could be associated with rainfall or earthquakes. At both stations, radon levels detected have variation lower than the radon standard deviation σ, in contrast to the Scholari data (Fig. 3).

Following the results after 2000, an inverse behaviour of radon level with temperature was observed at Gerakarou, while no influence appears with pressure, indicating that diffusion through soil porosity could be the dominant mechanism of soil gas transfer in a common fracture underground [16,17,18]. On the contrary, Stivos and Scholari's results show an inverse behaviour between radon and pressure, while temperature does not affect the radon level in soil gas, suggesting that advection through fractures and cracks in the underground is the primary process [15, 21]. Interestingly, at Sholari station, after the earthquakes occurred in the Autumn of 2001, the pressure and radon followed inverse behaviour, with pressure appearing similar to the Stivos station, where it was always detected. A possible origin of the lower variability of soil Rn concentration observed at Stivos station could be that the station is located at the junction of the two faulting directions presented in Fig. 1. In contrast, Gerakarou and Scholari stations are located at each one of the faulting directions, with the Scholari to be the closet station to the epicentre of the 1978 event.

However, it is more important and valuable to compare the results obtained from all three stations to evaluate the effects of earthquakes on soil radon concentration, but the lower than the standard deviations of radon variability in soil gas detected at Gerakarou and Stivos stations gives Scholari measurements the advantage for further time series analysis. The 15-min raw data collected are entered into a Python program and averaged daily to exclude the daily radon variation. After that, the daily data were averaged every 30 days (monthly). The program performed ANOVA analysis, and the difference between the daily and monthly results produced by ANOVA analysis was input into Mathematica to analyze the time series of radon, temperature and pressure curves, see Fig. 5. A fitting process was applied, defining the cosine curve as the standard, and the results were for \({\text{Radon}} = - 71.14 \cdot \cos \left( {1.06 \cdot t - 9.01} \right)\), for \({\text{Temperture}} = 0.015 \cdot \cos \left( { - 4.64 \cdot t + 1.01} \right)\) and for \({\text{Pressure}} = - 0.196 \cdot \cos \left( {0.99 \cdot t + 1.21} \right)\). The x-axis is expressed in days after the analysis started (14 February 2000). Significant variation in radon levels was observed during the Autumn of 2001 and 2002 when the pressure and temperature variations were the lowest (Fig. 5). Because radon levels at Sholari station seem to be influenced by the pressure, as illustrated in Fig. 3, the fit radon curve and pressure combined in a system with unknowns the time (t) and the radon level, which were estimated considering only the pressure known. The difference between the measured in the field and radon estimated values obtained by solving the system is presented in Fig. 5.

Fig. 5
figure 5

The difference in the time series between the daily and monthly averages according to ANOVA analysis for a radon, b temperature and c pressure. The x-axis is the day after the beginning of the monitoring (13 February 2000). The y-axis expresses the difference between the daily value measured and the monthly average calculated

The figures show the signals of radon drops caused by precipitation and seismic events at the peaks. The signals related to earthquakes, since the variations were 2–3 σ, the radon standard deviation of 5 kBq −3 in the Scholari station. A long and complex radon signal appeared with a drop observed after rains in July 2001, followed by an increase of 10 kBq.m−3 before an M = 3.9 occurred on 10 August 2001 (543 days from 14 Feb 2000), predicting 17 days before. Then, another drop of 10 kBq m−3 due to rainfall in late July, followed by a significant boost to 15 kBq.m−3, appears in a period of low rainfall. The variation was detected 23 days before an event swan M = 4.0 and 4.2 happened on 8 October 2001 (602 days from 14 Feb 2000). Another short radon variation has a clear peak of 15 kBq.m−3 detected 12 days back to ML = 3.8 occurring on 8 Mars 2002 (753 days from 14 Feb 2000) (Fig. 6).

Fig. 6
figure 6

The radon difference (Bq m−3) between the measured in the field and radon estimated values obtained by solving the system of radon curve and pressure combined with unknowns the time (t) and the radon level, which was estimated considering only the pressure known. The x-axis is the day after the beginning of the monitoring (13 February 2000)

The association of the seismic events with the radon signals is based on the relation proposed between the magnitude M and the epicentral distance R (km). The Dobrovolsky (1979) equation, log R = 0.43 M, transformed by Hauksson and Goddard [35] as M = 2.4 log R – 0.43 and by Morelli [36] as M ≥ (2.3 ± 0.2) log (R) − (0.4 ± 0.3) gives the maximum possible distance, R, where a radon station should be located from the epicentral to detect precursory radon signals. Earthquakes with magnitudes M3–4 can be detected by stations between 20 and 70 km away, as in the Lagadas basin during the study, while for M4.5 -5.5, the distance increased to 90–300 km.

The results of continuous radon monitoring support the data obtained using integrated measurements with SSNTD in the previous years (Table 1). Then, radon peaked 30 and 45 days after two earthquakes, M = 4.7 (August 26, 1983) and M = 5.2 (February 19, 1984), with epicentral distances of 15 and 150 km, respectively [26]. Also, radon signal registration 30 days before an earthquake ML = 4.5 occurred on December 12, 1999 [27], the signal remained high, while a secondary radon peak appeared 15 days before an earthquake M = 3.4 occurred on January 31, 2000, and after returning to the normal level. The epicentral distance was around 20–30 km. The predicting time estimated for events lower than M = 4 agrees with data detected in the Ioannina area, where predicting time ranged from 4 to 15 days before earthquakes with M = 4–3, and the radon variations were 3 times higher than the normal [33].

Table 1 Radon, soil gas signals, associated with earthquakes detected using integrated sampling and continuous sampling at the Stivos area

The results collected over two decades from the Stivos area allow comparison with equations proposed between precursor time (tp, days) and earthquake magnitude (M). The proposed equations suggest log tp, = A + B*M, where A and B are statistically determined parameters. Based on that, Rikitake 1976 calculated log tp, = − 1.83 + 0.7 * M [22], although [25] proposed log tp, = 0.32 * M [25]. The data obtained from continuous radon monitoring show a higher predicting time than the integrated measurements for the same magnitude (see Table 1). If only the continuous radon data were fitted, then the Rikitake 1976 equation seems more applicable. However, the Virk equation is superior when it fits all the data obtained from this area (Fig. 7) since an equation log tp, = 0.36 * M was fitted to the experimental data with R2 = 89%, almost identical to the Virk 1995 proposal. It has to be considered that the equations proposed above have been based on earthquakes with M > 5 but are still applicable to data associated with events with M < 5.

Fig. 7
figure 7

Data between prediction time (tp, days) and local magnitude (M) of earthquakes recorded in the Stivos area. Closed circles in blue and red were derived from integrated radon measurements and continuous radon monitoring, respectively. Dash curves in red and blue were obtained from Virk's [25] and Rikitake [23] relationships. The solid line represents the fitting to the data obtained

The specific equations could be referred to as every seismic fault. Still, each fault has a particular personality structure. Many explanations have been proposed regarding the earthquake occurrence, as the crack-avalanche model [1] proposes an underground weathering formation due to accelerating cracking in the region of a future earthquake focal area, so stress corrosion in the rock matrix may be associated with abnormal behaviour of radon level. A similar proposition was published based on the dilatancy diffusion model; the earthquake forces crack the rocks, increasing the total porosity and radon entries into the porous, giving a signal under steady conditions. At the same time, pressure decreases in the entire preparation zone due to increased cracking and radon driven to the surface [9].

Based on this theory, another proposed equation includes, among the others, the relative shift detected in radon signal level ΔR = (Rp − Rb)/Rb, where Rp is the radon level of the signal, and Rb is the average radon level. Then, considering the radon lifetime, the following equation gives the prediction time before earthquake magnitude ML = 2 * log (ΔR/α*T) − 15.3, where α is statistically estimated, including the radon decay constant [23]. The seismic events at the Stivos faulting, North Greece, associated with radon signals detected by integrated and continuous radon monitoring, as presented in Table 1, with radon shift from 2 up to 6 times, demonstrate that the value α is around 9.1 × 10–11 d−1. The equation results are illustrated in Fig. 8. This value approaching 1 × 10–10 d−1 is one order of magnitude higher than that estimated using data from the Himalayas and reported as 1.8 × 10–11 d−1 [23], maybe due to different geomorphological structure and seismological behaviour between the two areas, one sited in a strong energy density (J m−3) region while the other in a lower one.

Founded on the Ramola et al. [23] proposed equation (Fig. 8), low radon signals were detected long before the event corresponds to low magnitude (< 4). At the same time, high radon signals occur shortly before a catastrophic earthquake. The Stivos data shows 2 to 6 times the average radon registered 12- 45 days before M = 3.8–5.2 events occurred. If a radon signal > 12 times the average detected, catastrophic earthquake M > 8 shall appear in < 1–3 days or severe M = 6–8 in 5–1 weeks, respectively. To resolve that issue and handle extremely high radon signals, long-term radon monitoring for a decade or more regarding improving this method. Moreover, more stations should be employed to improve statistics since radon varies through space and time significantly.

Fig. 8
figure 8

The prediction time in days detected for earthquakes with magnitude M in the Stivos faulting associated with radon signal shift (ΔR, in times above normal), based on [23] equation [25]. The colours represent the earthquake magnitude (M) associated with the radon signals

Conclusions

A time series analysis was conducted at the continuous radon level data, temperature, and pressure to find clear radon signals from rainfall and earthquakes. The background radon level was 25 ± 5 kBq m−3 at the Gerakarou and Sholari stations. A higher level of 50 ± 15 kBq m−3 appeared at Stivos, at the junction of faulting directions in the area. Clear radon peaks appeared due to 50–100 mm rainfall during a few days. A radon anomaly associated with events M 3.8–4.2 was detected 12 up to 36 days before. The data was complemented with previous results obtained from the same area. The seismic events were associated with the radon signal following that the epicentral distance from the radon station was in the frame of the well-adopted equations previously published. These and previous data of the prediction time (tp day) associated with earthquake M following the equation log T = 0.36*M, almost identical to the [25] proposal. Since the radon signal level ΔR considered, the [23] equation M = 2 * log (ΔR/α*T) − 15.3 fits the data well with α = 9.1 × 10–11. The Stivos data shows 2 to 6 times more than average radon levels registered 12–45 days before ML = 3.8–5.2 events occurred. However, more experimental data must be associated with earthquakes to reach safe conclusions.