Introduction

Radon is a chemical element of the VIII group that exists in the gaseous state under normal temperature and pressure conditions: its solubility in water and organic solvents increases with decreasing temperature (230 cm\(^3\)/kg at 20 \(^\circ\)C and 1 atm)1. The most relevant property is its radioactive behaviour, as detectable in various isotopes within the decay series of \(^{238}U\) (\(^{222}Rn\)), \(^{235}U\) (\(^{219}Rn\)), and \(^{232}Th\) (\(^{220}Rn\))1,2,3. In the natural environment, radon can be associated with rock formations affected by the presence of uranium and thorium, from which it is mobilized according to the elastic rebound model1,2. After decay, the expulsion of the radon atom from the crystal lattice of the precursor radioelement (\(^{226}Ra\)) occurs in the opposite direction to the emission of the alpha particle, released in the meantime. In details, \(^4_2\alpha\) is positively charged and the displacement must be at least 10 nm at an energy of 100 keV to verify this phenomenon4. Only \(^{222}Rn\) has a sufficiently long half-life (approximately 3.86 days) to be significantly detectable in the natural environment: the mean free path of the atom due to rebound has been estimated at 0.1 \(\upmu\)m in air and 63 \(\upmu\)m in water5. Transport involves ascending diffusion through fractures, openings, or porosities in lithological formations, sometimes increased by convective processes in association with gases such as \(CO_2\), \(CH_4\), and \(N_2\)6. Radon can further transfer into the aqueous phase, although the low value of the solubility coefficient (0.25 at 20 \(^\circ\)C) reveals a clear tendency for atmospheric accumulation2,7. Nevertheless, referring to aquifer systems, the gas undergoes advective transport phenomena, being affected by low chemical reactivity as a noble element: therefore, it can be assumed that water contributes to removing radon from the lithosphere, subsequently releasing it into the atmosphere8,9. In literature, the study of radon mobilized by water circulation is generally aimed at evaluating the recharge dynamics of karst aquifers under varying hydrological conditions10,11,12,13. On the other hand, several works describe some atmospheric processes in connection with radon, influenced by hypogeal morphology and differences in terms of pressure or temperature between the internal and external environments14,15,16,17. Additionally, hypogeous streams in karst systems can exert an additional dragging effect on air masses17. Dead-end caves exhibit a series of common characteristics regarding the relationship between atmospheric circulation and radon concentrations. Air presents two seasonal flows depending on temperature gradients, resulting in different amounts of radon recorded by instruments (low in winter and high in summer). In the former case, the flow of cold and dense air to the cave floor confines the gas in the fractures of the rock mass, while in the latter, the increase in temperature gradients between internal and external environments allows its diffusion into the cave and the removal by convective circulation, entering at the ceiling and exiting at the floor in summer15,18,19,20,21. Lithospheric radon flows can also be related to seismic events1,2,4,6,22, in addition to indicators such as fluctuations in the electric or magnetic fields and the ratio between velocities of compressional (\(v_p\)) and shear (\(v_s\)) seismic waves23, or \(CO_2\) behaviour24. Numerous site-specific empirical equations correlate magnitude with the observation time of the precursor phenomenon25,26,27,28. The identification of an anomaly in the continuous radon signal compared to background noise is mainly determined through statistical analyses29,30,31,32. Attempts to highlight correlations have been conducted since the 1960s2: in the Alpine arc, the most interesting results were obtained in the Carnic Alps (Friuli, NE Italy), where correlations were revealed both in the short term for earthquakes with \(M>3.0\) and in the long term for \(M>5.0\)33,34. Similar conclusions were inferred in the Maritime Alps from karstic springs observation in relation to local earthquakes35. The main purpose of this research is to analyze the continuously acquired radon concentrations in air (starting from 2003) and in water (starting from 2010) inside the Bossea dead-end cave. These data are cross-correlated with discharge values referring to the main drain of the karst aquifer and a secondary seepage, in order to assess the connection between the radioelement concentrations and the hydrodynamic behaviour. Finally, After detrending the radon emissions to eliminate the incidence of the water flow, the residual values are analysed to assess the potential presence of anomalies linked to the local seismic activity.

Test site

Bossea dead-end cave is located in the Corsaglia valley (Piedmont, Italy) near the transition from the Western Ligurian Alps to the Maritime Alps (Fig. 1a). Access is located at 836 m. The cave extends for approximately 3 km up to 964 m of elevation corresponding to the terminal sumps36. This area is part of the Internal Ligurian Briançonnais domain, with Permian volcanic and volcanoclastic rocks covered by a partially detached Meso-Cenozoic sequence lying on a Variscan basement located at significant depths36,37. These formations have undergone metamorphism during the Alpine orogeny, with temperature conditions of around 350 \(^\circ\)C and a maximum pressure of 0.4 GPa, attributable to the greenschist facies (chlorite-albite paragenesis)36. Consequently, a bottom-up outcrop of meta-volcanic rocks, quartzites, argillites, and meta-carbonates (marbles) is observable (Fig. 1b); therefore the cave can be divided in two different sectors (Fig. 1c):

  • The upper canyon segment, where the main underground drain of the karst aquifer (named Mora creek) runs through a canyon carved into Mesozoic marbles with several relict phreatic conducts36;

  • The lower halls and massive collapse segment, extending from Ernestina waterfall to the semi-abandoned tunnel used for access and reactivated only in cases of exceptional flooding38, where the contact between the Mesozoic meta-carbonates and the Permian meta-volcanics can be observed37;

Fig. 1
figure 1

Morphogeological asset of Bossea cave, with indication of measurement sites for radon and water flows. (a) Location of the investigated site. Basemap credits: Google Satellite WMS - QuickMapServices plugin edited in QGIS v. 3.34.11 (https://nextgis.com/blog/quickmapservices/). (b) Geological map of the in Bossea area. 1: Metarhyolites, metadacites, porphyroids (Permian); 2: Argillaceous schists and quartzites (Lower Triassic); 3: Marble limestones, metadolomites, calcareous schists, and micaschists (Middle Triassic); 4: Varicolored pelites. Massive, laminated, and brecciated dolomites with cavities from the dissolution of evaporitic minerals. Dark micritic limestones (Middle-Upper Triassic); 5: Limestones and marbles (Middle-Upper Jurassic); 6: Locally marbleized calcareous schists (Upper Cretaceous); Basemap credits: https://webgis.arpa.piemonte.it/ags/rest/services/geologia/Geo_Piemonte_250k/FeatureServer. (c) Positions of radon measuring instruments in air and water within Bossea cave, as well as those acquiring water flow rates and geochemical parameters data (electrical conductivity and water temperature). Modified from vertical long profile East-West defined by Club Alpino Italiano, S.O. Bossea.

Deformation and brittle fracture phenomena can be referred to a transpressive stress state during the Oligocene, causing isoclinal folding with axes oriented NW-SE in the western sector of the karst system and WNW-ESE in the eastern one36. The lateral contact between heterogeneous formations is observable in both some secondary reverse faults oriented parallel to the fold axis (contact between carbonates and quartzites with offsets of meters) and at the two main WNW-ESE strike-slip faults delimiting Bossea aquifer to the north (Prel Line) and south (Fontane Line, involving marbles and meta-volcanics). Furthermore, a NE–SW fault with mixed normal and left-lateral strike-slip behavior outlines the detachment plane in the lower sector36. Chemical dissolution process mainly shaped the upper marbles of the cave. On the contrary, erosion process almost exclusively affected the antiform of insoluble cataclastic meta-volcanics, as water infiltrating into the fractures of the carbonate block reached the detachment surface36,37,39. The deepening of the cavity caused gravitational instability in the meta-carbonate duplexes at the roof, involving them in numerous collapse processes that shaped the large halls of the lower section; subsequently, the collapsed blocks were subjected to a dual process of erosion and dissolution, predominantly by the action of the Mora creek36.

The surface morphology is characterized by steep SE facing slopes, with bare rock outcrops and mostly thin soil cover, as well as broadleaf forest (laryx) alternating with mountain meadows39. Alpine climate is influenced by the mitigating action of the Mediterranean Sea: rainfalls mainly occur in autumn and spring as well as snowfall in the winter period40. The recharge area of the aquifer covers a surface of 6.5 km\(^2\): altitudes vary between 800 and 1800 m40,41,42. This complexity is reflected in different supplies of the aquifer:

  • Autogenic recharge (primary inputs), involving rainfalls and snowfalls which infiltrates through Middle Triassic–Cretaceous limestones and dolomitic limestones, covered by eluvio-colluvial deposits and confined by impermeable quartzites and metavolcanics37,41;

  • Allogenic recharge (secondary inputs), involving impermeable pre-Triassic formation which can feed the aquifer through sinkholes located along detachment surfaces37,39;

Bossea aquifer is thus characterized by high permeability: the Mora creek receives contributions from both the carbonate-dolomitic massif and the metamorphic basement rocks41,42. Additionally, meta-volcanics can also supply a series of secondary seepages known as Polle and Stillicidi due to their secondary porosity (Fig. 1c). The former are distinguished by open and karstified fractures at the contact between marbles and the underlying basement; flows are below 1 l/h, as detectable in sites such as Polla delle Anatre, Polla dell’Orso, and Polletta. The latter are located on the cave’s vault and appear as fractures covered by abundant calcite deposits, with flows below 0.008 l/s (Milano, Torre, Sacrestia)41,42,43.

Radiometric analyses (\(\gamma\)-spectrometry) of samples from Bossea cave conducted in 2008 revealed higher amounts of radioelements in meta-volcanics compared to limestones and quartzites, with intermediate behaviour in porphyritic schists. Hence, It can be assumed that these concentrations primarily occurred at the interface between the carbonate complex and the underlying impermeable basement44,45,46, followed by radioactive decay generating several isotopic species. These results set the premise for investigating the relationship between the Mora creek and Polla delle Anatre seepage flow rates with the radon values, which are monitored through various measurement points located inside the cave (see “Methods” section).

Results and discussion

Radon and water

A temporal series between September 2020 and May 2021 describing the Mora creek behaviour is showed in Fig. 2a: the main peaks in water radon concentration can be detected after important flood events, such as the exceptional one occurred on 3 October 2020 (684 l/s), followed by a maximum on 4 October of 1633 Bq/m\(^3\). On 4 April 2021, 1375 Bq/m\(^3\) are detected after the measurement of 332 l/s the day before (3 April). While the former represents the response to intense rainfalls into the karst network, this latter flood event is attributable to snowmelt occurring in spring due to the increase of external temperatures (Fig. 2b). A delay of 46 h for radon is identifiable with AlphaGuard instrumentation by applying cross-correlation function (Fig. 3a,b and Table 1): flows are regulated by neo-infiltration contributions, which apply the hydraulic pressure necessary to mobilize long-resident volumes (piston flow model37,40), and thus with higher radon concentrations. The same aspects are also detectable by monitoring of other geochemical parameters38,39,41,42 (Fig. 4a). However, cross-correlation yields to an average estimation of the delay of the radon peaks, considering that the lag slightly varies according to the saturation conditions of the karstic network when referring to a specific flood event.

Fig. 2
figure 2

(a) Radon data acquired by AlphaGuard instrument in water (Torrente Mora) and cave atmosphere from Sep 2020 to May 2021: flood maximum value is reached on 3 Oct 2020 (684 l/s) due to rainfalls (293 mm on 2 Oct 2020), while that on 3 Apr 2021 (332 l/s) is attributable to snowmelt. (c) AlphaGuard data related to Polla delle Anatre site from Jan 2015 to Jul 2015: water peaks are referable to heavy rainfalls which generated a marked piston flow hydrodynamic response. (e) Data acquired by Sun Nuclear instruments from Oct 2008 to Aug 2009. Two peaks in water signal are recorded on 5 Nov 2008 (488 l/s) and on 28 Apr 2009 (1079 l/s) due to rainfalls (101 mm on 2 Nov 2008 and 126 mm on 27 Apr 2009 respectively); snowmelt is observable in early April through small oscillations in radon signal at the Upper Laboratory. Daily rainfalls and external average temperatures for each event are available on ARPA Piemonte database https://www.geoportale.piemonte.it/ and plotted in (b, d, f). AG—AlphaGuard; SN—Sun Nuclear; ATM—atmosphere; AQ—water; ST—Sala del Tempio; UL—Upper Laboratory; TM—Mora creek; PA—Polla delle Anatre; Ext T.—average external daily temperatures.

Fig. 3
figure 3

(a) Cross-correlation (\(\rho _{xy}\)) results for a maximum lag of 400 days. (b) Zoom on lags lower than 120 hours referring to AlphaGuard instruments. (c) Zoom on lags lower than 10 days referring to Sun Nuclear instruments. The symbol (*) refers to the maximum value of \(\rho _{xy}\) in each diagram: numerical values are listed in Table 1. AG—AlphaGuard; SN—Sun Nuclear; ATM—atmosphere; AQ—water; ST—Sala del Tempio; UL—Upper Laboratory; TM—Mora creek.

Fig. 4
figure 4

Electrical conductivity and in-water temperature values measured in the (a) Mora creek and (b) Polla delle Anatre seepage. Comparing these indicators with flow rates and external meteoclimatic parameters, it is possible to observe the piston flow hydrodynamic response of the karst aquifer. It is noteworthy to underline the specificity of the Mora creek, as an opposite water temperature trend compared to the electrical conductivity can be observed in case of flooding. This aspect is justified by the high altitude of the aquifer recharge area feeding the main drain, while the flow path of the Polla directly develops in contact with the meta-volcanics, generating a standard piston flow response (electrical conductivity and water temperature variations are consistent). EC—electrical conductivity; Water T—temperature in water; Ext. T.—average external daily temperatures; TM—Mora creek; PA—Polla delle Anatre seepage.

Table 1 Cross-correlation values referring to water or atmospheric radon and the Mora creek flows. In the former case, \(\rho _{xy}\) is the least accurate (0.74) due to a not negligible singularity affecting each flood event, which is strongly influenced by saturation conditions of karstic network and reflects on radon advective transport. Values of \(\rho _{xy}\) in atmosphere result higher, even though the presence of more or less significant data gaps and the temporal extension of the observation windows may fluctuate in the range from  0.82 to 0.95. Lower lags are detectable in the upper segment of the cave (between 1 and 2 days) if compared to the lower halls, affected by delays greater than 48 h.

In the Polla delle Anatre site, excess of data gaps does not allow for the application of cross-correlation function. This occurred for a variety of technical reasons, such as intrumental drifts, malfunctioning and lack of power supply. This site suffered from these problems more than other monitoring points mainly due to the difficulty in the quantification of radon present in the water at exteremely low flow rates (secondary seepage). Therefore, delay is manually calculated by comparing peaks recorded between January and June 2015, where uninterrupted data prevail (Fig. 2c). For example, a radon minimum of 1328 Bq/m\(^3\) is detected on 27 March after a flow peak of 1.3 l/s occurred on the same day with a lag of 6 hours. In this case, flood events are responsible for a decrease in radon concentrations, as the supply circuit of this secondary seepage develops in direct contact with the impermeable radioelement-rich basement37,39,42. Increase in flow rate associated to significant flood events coincides with the contribution of water from meta-carbonates, thereby diluting the radon emissions. Furthermore, the piston flow response of Polla delle Anatre results in a shorter lag between the positive flow peak and the radon minimum, suggesting the transmission of a larger pressure wave within fractures if compared to the Mora creek, likely generated by rainfalls and not by snowmelt, which do not affect this water flow (Figs. 2d and 4b). Full temporal series of radon data in water (where available) are reported in Supplementary Material A.

Radon and atmosphere

Radon in atmosphere exhibits a trend similar to the one recorded in water: delays vary according to the location of instruments inside the cave. In the Upper Laboratory concentrations are the highest overall: for example, as shown in Fig. 2a, a peak of 1314 Bq/m\(^3\) is recorded on 4 April 2021 (the day after the flood event), whereas in Sala del Tempio this maximum is reached only on 5 April (1053 Bq/m\(^3\)). Cross-correlation reveals a lag of 45 h in the former site and 66 h in the latter (Fig. 3a,b and Table 1). In the upper segment there is a lower atmospheric volume if compared to the large halls at the contact with the meta-volcanic rocks, where the radon diffusion process requires additional time to reach equilibrium with that transported by the aqueous phase. Indeed, in the upper canyon radon concentrations are the result of a dual effect in the atmosphere resulting from the mobilization of water volumes (piston flow) along with a significant flux of air mobilized from the innermost sectors of the cave (even preceding radon transported by advection in water). Once the lower segment is reached (from Ernestina Waterfall downwards), this second contribution becomes irrelevant, allowing only turbulent motion in water to mediate the release of radon into the atmosphere. This aspect is proved by air velocity measurements performed in lower halls, generally ranging from 0.06 to 0.45 m/s except near waterfalls where they increase up to 5–6 m/s46,47. Analogous lags can be identified by Sun Nuclear instruments (Fig. 3c), although acquired with lower frequency (Table 3); in this case, small atmospheric radon oscillations due to snowmelt are clearly detectable in early April (Fig. 2e,f), both in the upper segment (around 2400 Bq/m\(^3\)) and in the lower one (around 2000 Bq/m\(^3\)).

Unlike similar dead-end caves, these results demonstrate the absence of seasonal effects on regulating atmospheric hypogeal radon concentrations by convective cells. Bossea complex is characterized by a convective double-cell model, where only the outermost one is affected by seasonal changes48. The innermost one exhibits a constant behaviour throughout the year: the Mora creek is assumed to be the only factor regulating air circulation and consequently radon mobilization. Finally, in the Porfiroidi site (Fig. 1c), it should be possible to intercept radon directly diffusing from meta-vulcanites, but the high fragmentation of the dataset does not allow for detailed evaluations: it is only possible to notice a gas flow with an order of magnitude higher than the other measuring stations (10\(^4\) vs 10\(^3\) Bq/m\(^3\)). These data are fully reported in Supplementary Material A.

Radon and earthquakes; further hydrodynamic details

Incidence of earthquakes on radon peaks is determined after removing the effect of hydrodynamics on radon concentration and setting a critical threshold of 2\(\sigma\) to identify potential emissions related to local seismicity. In Fig. 5a a time window is set referring to two earthquakes occurred on 15 October 2013 02:47 (UTC +1) and 31 October 2013 22:09 (UTC +1). Residuals shown in Fig. 5b are determined from atmospheric (Sala del Tempio) and water (Mora creek) radon. Statistical parameters (Table 2) suggest poor correlations with seismic events: absence of significant peaks beyond the conventional anomaly threshold of \(+2 \sigma\) is shown; however, these results allow for further detailing of the aquifer behaviour. Considering a seasonal observation window, residual mean value (10\(^{-8}\) in water and 10\(^{-6}\) in air) suggests that in Mora creek positive values can be found under low-flow conditions. Normalization by imposing unitary areas underlying the curves implies that physical variables with a more homogeneous distribution over time (radon) occasionally assume lower values compared to a signal with significant peaks localized at limited instances (water flow). The amount of gas transported in water primarily depends on the solubility equilibrium typical of a chemically non-reactive substance, which remains almost constant due to the minimal fluctuation in the main drain temperature under low-flow conditions. In case of a flood event, piston flow effects become predominant: there is an increase in residuals in the first 48 hours, suggesting an increase in volumes of water mobilized rich in radon, followed by a decrease to significantly negative residual values in the following 24 h (exceeding \(-2\sigma\) threshold). The explanation for this second phase concerns with the arrival of water with low residence time in the phreatic conduits, hence with minimal radon concentrations. Obviously, during the recession period of the flood curve (variable duration up to 10 days), positive residuals typical of low-flow conditions are restored. All considerations are valid also for the Polla delle Anatre site, even though the higher content of radon in the aqueous phase tends to attenuate the negative minima just below the mean value. An example related to this site is shown in Fig. 5c using atmospheric radon monitored in Sala del Tempio: residuals are all below the \(2\sigma\) threshold (Fig. 5d). Same results are obtained when analyzing earthquakes with small magnitudes but reduced epicentral distances (Fig. 5e,f). Complete diagrams of residuals for earthquakes selected according to criteria exposed in “Methods” section are reported in Supplementary Material B.

Fig. 5
figure 5

Residuals obtained by eliminating hydrodynamic trends on radon signals in water and atmosphere. Mean values \(\mu\), standard deviations \(\sigma\) and \(2\sigma\) thresholds are shown in each time window. (a) Radon monitored in Mora creek and in the atmosphere of Sala del Tempio and its residuals (b); (c) Radon in Polla delle Anatre and atmosphere of Sala del Tempio with its residuals (d); (e) Radon in Mora creek and Sala del Tempio with residuals (f) considering earthquakes characterized by low magnitude but reduced epicentral distances from Bossea cave. For details about seismic events and values of statistical parameters see Table 2.

Table 2 Statiscal parameters determined by applying the residual method on radon concentrations acquired in water and atmosphere.

Conclusions

Radon constitutes an additional tool beyond the common geochemical parameters used to describe the hydrodynamics of a karst aquifer. In Bossea Cave, the monitoring results confirm not only the piston flow model that characterizes the water reservoir fed by autogenous and allogenic contribution. In addition, two distinct flow paths are distinguished between the main creek (Torrente Mora) and the various secondary inflows (Polla delle Anatre). Indeed, in the first case, an increase in radon concentrations in the aqueous phase is observed subsequently to the main flood events (with a lag of 46 h), due to the mobilization of water volumes with large residence time to ensure solubility equilibrium between water and radon diffusing in fractures from meta-volcanic rocks. On the contrary, in the secondary inflows, the rates remain constant except for exceptional flood events, where the influx of infiltrating water volumes from the meta-carbonate karst network results in dilution with the resident waters rich in radon (lag of 6 h), as this circuit is in close connection with the meta-volcanic rocks of the impermeable basement. However, further investigations are necessary in order to quantitatively validate the continuous measurements of radon in water.

In terms of atmospheric radon concentrations, there is a clear predominance of effects related to the features of underground hydrodynamics over seasonal fluctuations referable to the external environment. In the upper canyon of the cave, carved in the meta-carbonates, an increase in atmospheric radon is mainly detected 45 h after the beginning of Torrente Mora flood event, due to the piston action of significant air volumes from the innermost sectors of the cave, compounded by the increase driven by the transfer of what is transported in the aqueous phase in turbulent motion within the first 48–72 h. This second effect exclusively prevails throughout the lower halls: average lag is about 66 h.

The application of the residual method in order to eliminate the influence of hydrodynamics on radon trends does not reveal any anomaly to be significantly attributable to the mobilization of this radioactive gas in the lithosphere yielded by variations in local stress or deformation state induced by an earthquake. Indeed, the exceeding of \(-2 \sigma\) threshold is due to the unbalance between the volumes of radon and water which are mobilized during floods.

The results of this study would be of interest for applications that may be related to human exposure to radon, specifically in the matter of occupational safety. The reason lies in the fact that professional workers (e.g. speleological guides or laboratory technicians) are exposed to doses for a prolonged time, resulting in disease of the respiratory system. Furthermore, in order to properly assess the relationship between radon emissions and the weak local seismicity of the area, further investigations should be carried out using instruments with high-frequency data resolution. The replication of the proposed methodology for radon emission monitoring in different matrices and data analysis is desirable to point out the incidence of the site specific effects in different karst systems. Finally, the improved knowledge of the structure of the aquifer structure provides a better comprehension of how the karst system responds to changing hydrological conditions. This aspect has applications in environmental monitoring and in assessing the availability of water resources, including potential supply for civil use.

Methods

All graphs and analytical evaluations are performed by using MATLAB (version R2022a, https://it.mathworks.com/products/matlab.html) and therein functions where indicated. Radon measurements are available as spreadsheets (.xlsx files). Maps shown in Fig. 1 are created by using QGIS (version 3.22.13, https://www.qgis.org/it/site/).

Atmospheric radon sampling has been conducting in Bossea Cave since 2003, involving the use of a passive Electret Ion Chambers46,47 (Fig. 1c). The working principle is based on a Teflon electrode that serves both as the source of the electrostatic field and as a sensor. Once diffused inside, radon is ionized and deposits on the surface of the electrode: concentration of the gas in the measured environment is determined in terms of the accumulated charge2,5. Timeline referring to positions of the instruments in the cave is reported in Table 3: devices produced by Saphymo GmbH, commercially known as AlphaGuard, provide a higher sampling frequency if compared to Sun Nuclear ones. Continuous measurement of radon concentrations in water (Table 4) involved the installation of experimental devices by Saphymo GmbH in 2010, consisting of a gas-permeable spiral tube screwed onto a rigid support, which is inserted into the water flow to be analyzed; a connected pumping system ensures a constant flow of radon to the AlphaGuard monitor. Response of radon emissions to the water flow are analyzed disregarding the absolute values that might be discrepant between different recording instruments adopted in the present and past monitoring periods46,47.

Table 3 Working period of instruments for continuous measurement of radon in cave atmosphere.
Table 4 Working period of instruments (located in the Main Laboratory) for continuous measurement of radon in water.

Radon atmospheric data are cross-correlated with measurements of water levels in the Mora creek, obtained through a rectangular reinforced concrete weir and associated to a digital data logger with a sampling of four hourly measurements since 2015 (acquiring also values of electrical conductivity and water temperature)37,39. The hydrological level (L) thus obtained needs to be converted into a flow rate Q using the classical hydraulic formula of the rectangular weir (1), given its width (b = 1.20 m) and the discharge coefficient (\(\mu =0.35\))48,49.

$$\begin{aligned} Q=\mu \cdot b \cdot L \cdot \sqrt{2gL} \end{aligned}$$
(1)

Flow rate Q of Polla delle Anatre is determined by using a triangular weir with digital data logger acquiring the level L every 15 min starting from 201537,39; the calculation is performed by using a formula defined specifically for this secondary flow (Eq. 2) with manual calibration48.

$$\begin{aligned} Q=41,392 \cdot (L-0,022)^{1,4273144} \end{aligned}$$
(2)

Seismic data are acquired from the INGV database updated in real-time (https://terremoti.ingv.it/), reporting earthquakes occurring since January 1, 1985 regardless of the magnitude value. Only epicenters located at a radial distance not exceeding 50 km and with \(M \ge 3.0\) are considered, coherently with studies conducted in close and analogous areas in terms of geological context35. For completeness, seismic events with \(M < 3.0\) but short epicentral distance (lower than 10 km) from Bossea cave are separately examined. These data are fully described in Supplementary Material B.

Residual method procedure aims to separate hydrodynamic effects from the radon signals acquired on site. Due to the low frequency of data acquisition and the presence of numerous temporal gaps, only measurements taken by AlphaGuard instruments were considered for this purpose, while Sun Nuclear ones were discarded. Residual values are defined by following these steps:

  1. 1.

    Selection of earthquakes that respect the epicentral distance and magnitude requirements;

  2. 2.

    Definition of the observation time window, selected on a seasonal basis in order to highlight the appropriate deviations of radon signal from the long-term defined average values;

  3. 3.

    Calculation of temporal lag between peaks of the water flow and radon concentration curves (using cross-correlation), followed by the subsequent translation of the former signal to align the peak with that of the gas (bringing it into the temporal reference system of the latter);

  4. 4.

    Normalization of radon and flow signals by imposing a unit value for the areas under the curves \(\int _{t=t_0}^{t=t_1}{f(x) \cdot dx =1}\). This criterion was chosen to minimize effects of signal gaps;

    $$\begin{aligned} f_{norm}(x)=\frac{f(x)}{\int _{t=t_0}^{t=t_1}{f(x) \cdot dx}} \end{aligned}$$
    (3)
  5. 5.

    Computation of the residuals by subtracting the flow signal from that of radon, assuming that they should not depend on the hydrodynamics of the karst system;

    $$\begin{aligned} residuals=Rn-Q \end{aligned}$$
    (4)
  6. 6.

    Analysis and interpretation of residuals, by determination of the mean value and standard deviation;

The general concept of anomaly in radon flux can be defined as the identification of a certain time interval in the continuously acquired signal during which gas concentrations exceed a given threshold. This trend can vary in duration, ranging from a few hours to several years, and can be either positive or negative relative to the background values29. Currently, there is no rigorous definition of the anomaly value: in this study a well-established practice is taken, referring to twice the standard deviation added to (or subtracted from, in the case of the lower limit) the mean value of the analyzed time series (\(2\sigma\) method2,30,31,32).