Evaluation of tectonically enhanced radon in fault zones by quantification of the radon activity index

This work highlights the importance of the Geogenic Radon Potential (GRP) component originated by degassing processes in fault zones. This Tectonically Enhanced Radon (TER) can increase radon concentration in soil gas and the inflow of radon in the buildings (Indoor Radon Concentrations, IRC). Although tectonically related radon enhancement is known in areas characterised by active faults, few studies have investigated radon migration processes in non-active fault zones. The Pusteria Valley (Bolzano, north-eastern Italy) represents an ideal geological setting to study the role of a non-seismic fault system in enhancing the geogenic radon. Here, most of the municipalities are characterised by high IRC. We performed soil gas surveys in three of these municipalities located along a wide section of the non-seismic Pusteria fault system characterised by a dense network of faults and fractures. Results highlight the presence of high Rn concentrations (up to 800 kBq·m−3) with anisotropic spatial patterns oriented along the main strike of the fault system. We calculated a Radon Activity Index (RAI) along north–south profiles across the Pusteria fault system and found that TER is linked to high fault geochemical activities. This evidence confirms that TER constitutes a significant component of GRP also along non-seismic faults.

Radon transport over distance in fault zones may locally represent a significant additional contribution to the Rn background originated by radionuclide decay in the source rocks; this quantity can be defined as Tectonically Enhanced Radon (TER). TER measures the "strength" of the tectonic factor that can potentially increase the radon originated from the local lithology (and, as a consequence, can affect the potential increase of the IRC), due to: -the action of carrier gases (e.g. CO 2 , CH 4 ) migrating from deep sources by advection along faults 7 .
Numerous studies report examples of radon migration along faults that may be ascribed to the TER concept but specific studies aimed at quantifying this geogenic component of radon are missing. However, most of these studies deal with seismic faults, in contrast to the aseismic fault system discussed in our work. Both seismic and aseismic faults can still be expected to have certain features in common; while aseismic faults can be supposed to lack the temporal dimension of seismic faults which lead to temporally variable gas transport and exhalation. This effect is used for diagnosis and forecast in volcanology and seismic science, but not further discussed here. For example, Ciotoli et al. 7 investigated radon distribution in the Fucino Plain (central Italy), a tectonically active intermontane basin bordered and crossed by a complex network of deeply buried and/or shallow faults and fractures characterized by high seismic activity (e.g., magnitude 7.0, Avezzano earthquake of 13 January 1915 15 ). In this area, radon anomalies up to 5 times the background soil production (25 kBq·m −3 ), occur along the main faults of the basin even when buried under thick sedimentary covers (up to 900 m). Fault-related anisotropy affects radon distribution at the surface, and provides a clear correlation between the shape and orientation of radon anomalies and the geometry of the recognised damage zones in the area. Seminsky et al. 11 studied the variations of radon activity in the crustal fault zone of the Baikal-Mongolian seismic belt and revealed that the radon anomalies occurring along the investigated profiles crossing the Khustai fault are five times higher than the background value (3.9 kBq·m −3 ). Wang et al. 16 found radon concentrations in soil gas eight times higher than the background value (4.7 kBq·m −3 ) along the active Tangshan fault (Northern China). In particular, radon enhanced values in the area showed close relation with the seismic activity along fault zones. Active segments of faults and associated damage zones along which larger earthquakes nucleated may act as preferential paths for enhanced radon migration. These studies suggest that the anisotropic spatial distribution of radon anomalies along fault zones is strictly related to their orientation, geometrical complexity (i.e. single fault density vs complex fault system), seismic activity and the relation between width of the core zone and extension of the damage zone (e.g., presence of transverse faults, i.e. relay ramp 17 ) (Fig. 1). with permeable core and damage zone, TER shallow signal is represented by a single-peak anomaly; (B) single fault with non-permeable core and wide and permeable damage zone, TER signal may show a double-peak anomaly; (C) complex fault zone with permeable core and damage zone, TER signal occurs as a single-peak; (D) complex fault zone with non-permeable core, and wide and permeable damage zone; TER is represented by a double-peak signal. In this work, we focus on the analysis of the distribution and the magnitude of radon, thoron and CO 2 soil gas anomalies in the Pusteria/Pustertal Valley (north-eastern Italy, Fig. 2). This area is characterised by the Periadriatic fault system which is the aseismic tectonic boundary between the Austroalpine orogenic wedge and the Southalpine indenter in the Eastern Alps. In the study area, this non-seismic fault system is composed of three main E-W trending faults: (i) the Pusteria/Pustertal fault (PF), a sub vertical fault system that accommodated a dextral transpressive kinematics during the Alpine orogenesis since late Oligocene 18-21 ; (ii) the Kalkstein-Vallarga/Weitental fault (KV), a minor fault still pertaining to the same Periadriatic fault system 22 ; (iii) the Deffereggen-Anterselva/Antholz-Valles/Vals line (DAV), a mylonitic shear zone marking the southern boundary of the Alpine (Paleogenic) metamorphic overprint within the Austroalpine basement 23 . In the study area, the outcropping crystalline basement is composed of orthogneiss and paragneiss (Austroalpine domain) to the north of the PF, and phyllites intruded by granites (Southalpine domain) to the south of the PF (Fig. 2).
The main objectives of our work are: (i) to evaluate the potential degassing processes along a non-seismically active fault (i.e., the Pusteria fault system), (ii) to quantify TER in terms of a Radon Activity Index (RAI) and iii) to compare RAI magnitude with those calculated for seismically active fault systems from the literature. The obtained results highlight that TER may constitute an important additional component to GRP also along aseismic fault systems, thus potentially enhancing influx of radon into buildings.

Results
222 Rn and thoron ( 220 Rn) activities, CO 2 and O 2 concentrations were measured directly in the field at 278 sites according to a regular grid and along three N-S profiles crossing the aseismic Pusteria fault system, over an area of about 60 km 2 . Table 1 shows descriptive statistics (i.e., main numerical indexes) of all the studied gases. The significance of Kolmogorov-Smirnov test (p < 0.01) confirms the non-normal distribution of all the variables (measured on the grid and along the three profiles), especially those of 222 Rn, 220 Rn and CO 2 which showed a skewed distribution (see Supplementary Fig. S1a).
Descriptive statistics on the total dataset (A) highlights that 222 Rn ranges between 1.54 and 815.96 kBq·m −3 . We considered the similarity between median (36.52 kBq·m −3 ) and geometric mean (GM) (35.16 kBq·m −3 ) as a representative indicator of the approximately log-normal distribution of the variables, since median is not influenced by the presence of outliers and GM is the mean value of the log-transformed data. 220 Rn values range The normal probability plots (see supplementary materials, Fig. S1b) describe both the background values and the anomaly thresholds of the different gas species: 50 kBq·m −3 for 222 Rn, 15 kBq·m −3 for 220 Rn and 4% v/v for CO 2. We applied the Akerblom formula 26 (see section on Radionuclide content in Methods) to calculate the radon at equilibrium with the activity concentration of the 226 Ra (its direct parent radionuclide) in the main outcropping lithologies. The formula provides 222 Rn mean values of 83.5 kBq·m −3 for the Austroalpine gneiss, 68.3 kBq·m −3 for the Permian granite and 27.7 kBq·m −3 for the Southalpine phyllite, respectively (see Supplementary Table S1). The mean value of 59.8 kBq·m −3 agrees with that obtained by the NPP (Normal Probability Plots; see Supplementary Fig. S1b).
Variogram analysis has been performed on 222 Rn, 220 Rn and CO 2 data which were log-transformed to minimize the effect of outliers [27][28][29] . The experimental variograms were calculated along four directions to highlight the presence of anisotropies (e.g., fault related gas anomalies) in the spatial distribution of the data. In general, all gases showed a clear anisotropic behaviour along the E-W direction. The experimental variograms calculated along the maximum axes of the anisotropy ellipse showed correlations consistent with the main directions of faults and fractures. In particular, the main anisotropy directions are about 90° for 222 Rn, 50° for 220 Rn and 95° for CO 2 (see Supplementary Fig. S2). Variogram models were used within the ordinary kriging algorithm to predict values at unsampled locations and construct gas distribution maps (see section on Statistical and Geospatial analyses in Methods).
According to the maximum anisotropy axis discussed above, 222 Rn anomalies well fit the E-W trend of the Periadriatic system from the DAV line, to the north, to the PF, to the south (Fig. 3a). This trend also follows the direction of the brittle fracture zone between the main faults, which includes conjugate minor faults and fractures roughly parallel to the main faults 30,31 . The highest anomalous values (> 100 kBq·m −3 and up to 800 kBq·m −3 ) occur in the eastern sector (municipality of Falzes/Pfalzen). This area shows the highest radon values up to 815.96 kBq· m −3 . In contrast, the plain area south of the PF displays the lowest radon values, below the background (< 50 kBq·m −3 ).
The lowest 220 Rn values (< 30 kBq·m −3 ) have been measured in the fracture zone between the DAV and PF, while the highest ones occur to the south of PF corresponding to the western sector (Terento/Terenten municipality) (Fig. 3b).  www.nature.com/scientificreports/ axis = 95°) (Fig. 3c). The area south of the PF shows a CO 2 anomaly greater than 4% in the municipality of Chienes/Kiens; in contrast, the easternmost sector of the study area only reports low CO 2 concentration values. In general, the most fractured sector between DAV and PF shows higher values of 222 Rn with respect to neighbouring zones. It is worth noting that 222 Rn and CO 2 anomalies coincide along the PF, whereas CO 2 displays high values also south of the PF.

Discussion
We performed morphological analysis from the DTM (2.5 × 2.5 m Digital Terrain Model) in order to reconstruct a structural model of the area (see section on Structural Model in Methods). The model highlights a series of minor faults and fractures located among the main faults in the area (PF, KV line and DAV). This structural network includes the potentially permeable pathways for gas advection.
As widely reported in the literature, radon generated by 238 U decay in rocks and soil can escape by diffusion from the surface and accumulate/disperse at shallow depth in soils, or migrate by advection upwards at greater distance from deeper sources along preferential pathways, such as faults and fractures, transported by gas carriers (CO 2 , CH 4 , etc.) 7,25,[32][33][34] . Meteorological parameters (e.g., pressure, humidity, temperature) can also affect the distribution of radon in shallow soil, but they cannot explain the sharp anisotropy of the radon distribution along the fault strike 35,36 .
In this study area, 222 Rn and CO 2 anomalies generally well fit with the main orientation of the Periadriatic fault system. However, despite both gases show the same anisotropy ratio (R = 2), as observed in the experimental variograms (see Supplementary Fig. S2), 222 Rn anomalies show higher spatial continuity compared to CO 2 . As a result, the spatial distribution of 222 Rn and CO 2 concentrations show different patterns characterised by elongated anomalies, and aligned spot anomalies, respectively, as also elsewhere observed by Ciotoli et al. 7 .
The extent and the elongation of radon anomalies spatially changes from a sharp (narrowly elongated anomalies) pattern in the western sector (from Terento/Terenten to Chienes/Kiens municipalities) to a more diffuse pattern in the eastern sector (Falzes/Pfalzen municipality), where the KV is displaced by secondary orthogonal NNE-SSW trending faults and width of the fractured area between KV and PF increases. It is worth noting that the distance (≈ 2 km) among these secondary faults is consistent with the range of the experimental variogram of the radon data measured along the major axis of the anisotropy ellipse. This parameter represents the distance at which Rn migration along the fault system shows a continuous behaviour, and thus indicates the extension of the permeable sectors along the PF system that can be dissected by minor faults according to a series of relay ramps that may interrupt the continuity of this tectonic system. Accordingly, the length of the minimum anisotropy axis (≈ 1000 m) well fit with the N-S extension of the elongated radon anomaly occurring from Terento/Terenten to Chienes/Kiens municipalities. Furthermore, though CO 2 anomalies show a spotty distribution pattern, they are also consistent with the orientation of the PF and confirms the role of CO 2 as carrier gas for 222 Rn from deeper sources. A scatterplot between radon and CO 2 values (see Supplementary Fig. S3), occurring in correspondence of the anomalous area along the Pusteria Fault, shows that there is a significant linear relationship between the two gas species (p value < 0.0001) with coefficients R = 0.723 and R 2 = 0.523. This result confirms that CO 2 acts as a carrier gas for radon along the main fault. In the southern sector (Chienes/Kiens municipality and S. Sigismondo/St. Sigmund village), low 222 Rn (below the background value of 50 kBq m −3 ) and high 220 Rn activities suggest the presence of high permeable soils in the shallow environment, considering the very short half-life of 220 Rn (55.6 s). According to this hypothesis, in this area CO 2 anomalies can be linked to shallower factors (e.g., soil humidity, organic matter) caused by the presence of the Rienza/Rienz river and farmed fields.
The eastern sector (Falzes/Pfalzen municipality) displays the highest 222 Rn activities (from 100 to 800 kBq·m −3 ) measured in the area; these values are comparable with maximum 222 Rn activity (order of magnitude 10 2 kBq·m −3 ) measured in different seismic intermontane basins in Italy 25 and in China 13 . In this sector, 222 Rn anomalies show a more diffuse pattern though they are also characterised by local E-W anisotropy due to the presence of minor faults and according to the direction of the main fault system. In contrast, CO 2 spatial distribution does not show anomalous concentrations. This might be due to an increase in soil ventilation (i.e., CO 2 dilution) caused by the higher elevation of this zone (> 1800 m.a.s.l.), the reduced soil thickness and rock weathering 36 (i.e., the diffused presence of slope deposits consisting of large loose rock blocks). However, all these factors do not seem to affect 222 Rn emission in terms of TER.
In contrast, in the western sector (Terento/Terenten municipality) 222 Rn anomalies show a more defined anisotropic pattern extending towards the central sector (Chienes/Kiens municipality) where the 222 Rn anomalies are mainly restricted between the PF and the KV. However, this area is characterised by a less pervasive large-scale fracturing and thus radon measurements show lower values than those measured in the eastern sector (mean 222 Rn activities of 75 kBq·m −3 and 130 kBq·m −3 , respectively).
Despite the calculated mean value of 222 Rn activity at equilibrium with 226 Ra is slightly higher than that observed in the NPP, the contribution of the lithology is not high enough to justify the magnitude of the anomalies measured in the eastern sector (Falzes/Pfalzen area) and along the whole Pusteria fault. Therefore, we can conclude that all values above the anomaly threshold of 50 kBq·m −3 are reasonably supported by an upward migration along the fracture zone associated to the PF that promotes tectonically enhanced radon at the surface thus contributing to the increase of the GRP in terms of TER.
Based on the above considerations, we can confirm that the study area shows different pattern of 222 Rn, CO 2 and 220 Rn spatial distribution of the gas anomalies at the surface, primarily related to the extension of the fracture zone and the length of the dislocation, and secondarily to local geological factors. These parameters mainly govern the pattern and the magnitude of the observed fault-linked Rn anomalies, thus suggesting a different geochemical activity (i.e., the gas-bearing property) of the fault segments that makes up the PF system 10 www.nature.com/scientificreports/ In order to evaluate the correlation between TER and geometry of the PF system, we firstly investigated the relationship between radon anomalies and fault/fracture density map obtained by Kernel density algorithm (see section on Structural Model in Methods) in order to identify TER zones. Then we tried to quantify the TER zones in terms of radon activity index (RAI) (see section on RAI Calculation in Methods) 11 . Figure 4 shows a summary sketch map of the areas characterised by fault-linked Rn anomalies (e.g., TER). The map was obtained by reclassifying and overlaying radon anomalies (above 50 kBq·m −3 ) and fault/fracture density map (> 1.5 fault/km 2 ). The figure also highlights the distribution of the main and minor faults reported in the structural model (see section on Structural Model in Methods).
In particular, the map shows 4 zones characterised by different radon sources and magnitude: • the zone of lithological background (white area), statistically and geochemically defined by radon values below 50 kBq·m −3 ; this zone (about 34 km 2 ) is characterised by Rn mean value of 20.5 kBq·m −3 ; • the zone of high fault/fracture density (> 1.5 fault/km 2 ) but with radon values below the background (pale grey area); this zone (about 14 km 2 ), where the presence of non-permeable fault segments is supposed, is characterised by Rn mean value of 22.3 kBq·m −3 • the zone with radon activity higher than the background (silver grey area); this area, apparently not correlated to zones with high fault density, mostly occurs in the fractured area recognised within the fault system. This zone (about 8 km 2 ) is characterised by high Rn mean value of 111.5 kBq·m −3 that suggests the presence of radon advection processes which change toward the surface in a more diffusive (even laterally) character; this may explain the more pervasive distribution of radon anomalies; • the TER zone (about 12 km 2 ) (dark grey area) is characterised by very high fault-linked radon anomalies with mean value of 154.0 kBq·m −3 mainly occurring in the area comprised between the KV and PF.
As TER represents the radon contribution from faulted areas, the quantitative evaluation of the contrast between the magnitude of the soil-radon anomalies (Rn Max ) near the fault and the lithological background (Rn BG ) was carried out by calculating the RAI for each fault-related radon concentration peak recognised along 21 N-S profiles (P4-P24) extrapolated from the estimated Rn grid map (Fig. 3a). The ratio between these two values can be considered as a function of the geochemical activity of the fault, as well as of the degree of fracturing in the damage zone. In general, high RAI values correspond to high Rn absolute concentration values in the soil-gas. They are located above damage fault zones, and have great practical interest for the assessment of radon hazard in populated areas. The RAI values obtained in the study area provide a quantification of the geochemical activity of the PF system within the recognised TER zones. RAI values were grouped according to the classification procedure described in the RAI calculation section of the methods, and displayed as classed post map (Fig. 4).
Calculated RAI values (see Supplementary Table S2) range from 1.08 to 13.86 and are grouped in four classes according to the following classification 11 scheme: low geochemical activity (RAI < 2), moderate geochemical activity (2 < RAI < 3), high geochemical activity (3 < RAI < 6) and very high geochemical activity (RAI > 6). The classed post map of the RAI values highlights that 25 fault-related radon concentration peaks out of 40 (about 62%) are localized in the fracture zone between the two main faults (PF and KV), most of the values falling in the moderate (9), high (5) and very high (5) geochemical activity groups; values below 2 only occur in the area of Chienes/Kiens municipality and to the south of the PF. The distribution of the RAI values in Fig. 4 describes a specific pattern which follows the fault strike direction and extends over the fracture zone which is 5 km wide to the east and 2 km wide to the west, where the two main faults PF and KV tend to converge.
High to very high RAI values mainly occur in the eastern sector (Falzes/Pfalzen municipality). In this area, the highest RAI values (mean value = 4.8) correspond to high radon concentration peaks aligned along the direction of minor faults in the middle of the TER area (dark grey zones in Fig. 4) between the PF and the DAV fault.
The western sector (Terento/Terenten municipality) displays medium to high RAI values (mean value = 2.9) and the associated radon peaks are strictly located within the TER area occurring between the PF and KV. In this sector, RAI values of the Rn peaks show a good spatial continuity as suggested by the variogram model, and the radon anomaly in this case is controlled by increased gas permeability of the damage zone along this fault segment of the PF.
The central sector (Chienes/Kiens municipality) displays low RAI values (mean value = 1.6). Although located within the TER zone (dark grey TER area in Fig. 4), this area is characterised by a patchy pattern of the radon anomalies with fault segments along the main structure of the PF with radon activity slightly above the background values (< 100 kBq·m −3 ). These fault segments, bordered by secondary NNE-SSW tending transfer faults with sinistral relay ramp geometry, are characterised by relatively unfractured rocks covered by weathered soil that lowers the gas permeability thus having a high impact on the intensity of radon anomalies near the fault zone. These conditions determine the patchy character of soil gas anomalies (including radon) at the surface (Fig. 1d), and confirm that in this sector of the valley the PF system, despite of the high fracture density (pale grey area in Fig. 4), gas permeability is highly variable along its entire extension 7,32 .
DAV mylonitic zone does not show any specific Rn peak pattern because it is characterised by a very low gas permeability due to its massive mylonitic microstructure. Furthermore, since this fault is located along the northern sector of the study area, boundary effects due to the low sample density makes the estimate less robust due to the lack of neighbouring points [37][38][39] . According to the distribution of the TER zones and their RAI values, the Pusteria fault system displays different distribution patterns of the anomalies and variable geochemical activity. In particular, in the western sector, where the different fault lines of the PF system tend to converge, the fault system behaves as a simple single fault (e.g., A and B in Fig. 1), whereas in the eastern sector the PF system acts as a complex fault system (e.g., C and D in Fig. 1). The distribution of the TER zones, explained in terms of RAI, confirms that they are characterised by high to very high radon values when associated with high fault/fracture density.
The relative RAI values obtained along different N-S profiles reflect the contrast between the radon anomalies and the background. According to the described data, RAI values can vary along fault strike within an order of magnitude depending on the type and fracture density of the fault zone. The individuation of contrasting anomalies and the calculation of the Radon Activity Index (i.e. RAI) is crucial to quantify TER in a fault zone and interpret the data in terms of geochemical activity. The obtained results are in line with the RAI calculated for some seismically active faults reported in the literature confirming that also non-active faults may increase the GRP (Table 2). We can establish that the collection of soil gas radon measurements for mapping the TER zones, as well as the calculation of RAI values in fault areas, are of practical interest for local authorities in the assessment of radon hazard, especially in urban areas. Table 2. Values of RAI obtained by the ratio between the maximum value measured in case study the area and the background value over different areas in the world.

Conclusions
The main goal of this study was to evaluate the potential degassing processes along an aseismic fault system, the Pusteria fault system, in north-eastern Italy. To accomplish this objective, we have defined a further geogenic radon component (Tectonically Enhanced Radon, TER) that represents the contribution of tectonics to the Geogenic Radon Potential (GRP) of an area. Then, we have quantified the TER component by calculating the Radon Activity Index (RAI) that defines the geochemical activity of a fault in terms of Rn emission. The obtained results highlight the following conclusions: -the fracture zone of the Pusteria fault system, including the main faults PF and KV, and other minor faults, plays a fundamental role for gas ( 222 Rn and CO 2 ) migration towards the surface thus providing TER as an additional contribution to the Rn diffusing from the source rocks (i.e., lithological background); both these components account for the GRP of an area; -the summary map of the TER zones (Fig. 4) confirms that most of the areas characterised by high radon values are also associated with high fault/fracture density; -radon activity index (RAI) calculated at Rn peak values along N-S profiles represents a good proxy to quantify the geochemical activity (e.g., gas-bearing properties) of a fault zone. We recognised that TER areas are characterised by medium to very high RAI values. However, although the central sector of the PF system shows high fault/fracture density, the low RAI values suggest the presence of less permeable segments along the fault system; -the comparison of the calculated RAI with those calculated from seismic areas with active faults, reported in the literature, highlights that also aseismic faults can provide conduits for Rn migration towards the surface with the same order of magnitude of seismic faults; -since the radon high activity generated from source rocks represents background areas of potential radon risk for inhabitants, the quantification and mapping of the TER areas are important to better evaluate the potential risk due to the increased radon availability to influx within buildings. -the knowledge of these two components of the GRP is fundamental for the mapping of susceptibility of the territory at different spatial scales and in different geological scenarios, and can help policy makers to plan monitoring activities and take mitigation actions reducing the damage for the society. The future developments of TER concept will have important implications on the identification of radon priority areas (RPA) as required by the European Directive 59/2013 3 .

Methods
Structural model. The geological map 45  Radionuclide content. A total of 14 rock samples of the main outcropping lithologies in the study area were collected as follow: gneiss (7 samples) from the Austroalpine unit; phyllite (3 samples) and Brixen granite (4 samples) from the Southalpine unit. Prior to radioactivity measurement, the rock samples were ground and mechanically sieved (< 4 mm) to determine the activity concentration of natural radionuclides ( 238 U, 232 Th, 40 K) using high-resolution gamma-ray spectrometer in a suitable configuration. Rock samples were analysed using two p-type coaxial Hyper Pure Germanium crystal detector (HPGe), a PROFILE (Ortec-Ametek Inc.) with an extended energy range (20-2000 keV) and a GEM model (Ortec-Ametek Inc.) with an energy range 80-2000 keV. These detectors have relative efficiency of 20% and 38%, and resolution (FWHM) at 1322.5 keV of 1.9 keV and 1.8 keV, respectively. Both systems were calibrated for energy and efficiency using liquid standard solutions (Eckert and Ziegler Multi-nuclide standard solution 7501) in a jar geometry (diameter = 56 mm; thickness = 10 mm). Spectra of the rock samples were acquired for one day to optimize peak analysis and subsequently processed and analysed using the Gamma Vision-32 software package (version 6.07, Ortec-Ametek). 226 Ra was determined at 186 keV correcting the peak area by the 235 U interference according to the method proposed by 47 , under the hypothesis of secular equilibrium between 226 Ra-238 U and natural 235 U/ 238 U isotopic ratio. 238 U and 232 Th were then determined using the emission of their radioactive daughters 226 Ra and 228 Ac. Data validation and quality control was carried out analysing a series of certified reference materials (CRMs) in the same geometry as the unknown samples and calibration standard. The CRMs used, in the solid phase, were: IAEA-412 Pacific Ocean Sediment, UTS-3 (Canmet) and Dh1-a (Canmet). Minimum detectable activity (CR MDA ) was determined using the Traditional ORTEC method (ORTEC, 2003) with a peak cut-off limit of 40%. Conversion from specific activity (Bq·kg −1 ) to bulk elemental weight fraction was obtained with the following conversion factors 48  The radon background production (as well as the geochemical anomaly threshold) has been calculated using the 222 Rn activity at equilibrium with parent radionuclides ( 226 Ra) in collected rock samples, using the Akerblom formula 26 : where C Rn and C Ra are radon in soil gas (Bq·L −1 ) and radium in soil (Bq·kg −1 ), respectively, ε is the emanation power coefficient (dimensionless), ρ is soil density (kg·L −1 ), and n is the effective porosity coefficient (dimensionless). The emanation power coefficient has been calculated as one minus the ratio between the mean activity concentration of radon decay products 214 Pb and 214 Bi with respect to the parent 226 Ra 49,50 . We have assumed that ρ = 2.7 g cm −3 and n = 0.4. Soil gas sampling. Soil gas surveys were performed in July 2021 to ensure stable meteorological conditions. A total of 278 soil gas samples have been collected in an area of about 60 km 2 according to a 500 m × 500 m grid and along 3 north-south profiles (P1, P2, P3) crossing the fault zone with a sampling step of about 100 m, in the municipalities of Terento/Terenten (western sector), Chienes/Kiens (central sector) and Falzes/Pfalzen (eastern sector).
We collected the soil gas samples using a 6.4 mm thick-walled stainless-steel probe pounded in the ground by a co-axial hammer at a depth of about 0.8-1 m to minimize as much as possible the influence of meteorological parameters that may affect the air exchange at the soil-atmosphere boundary. Soil gas measurements were conducted using a portable multi-gas analyser (Draeger X-am 7000, Drägerwerk AG&Co. KGaA) connected to the probe with a silicon tube for the simultaneous analyses of carbon dioxide (CO 2 , range 0-100%), methane (CH 4 , range 0-100% LEL), hydrogen (H 2 , range 0-600 ppm), hydrogen sulfide (H 2 S, range 0-1000 ppm) and oxygen (O 2 , range 0-21%). Radon ( 222 Rn) and thoron ( 220 Rn) measurements were conducted using RAD7 (Durridge Company, Inc) alpha detector (± 5% absolute accuracy, and a sensivity of 0.25 cpm/(pCi/L) 0.0067 cpm/(Bq·m −3 ). Each measurement of radon and thoron activity is performed with 5-min integration time, and is repeated until the difference between the last two measurements is at least below 5-10%. The final result was determined by taking the average of the last two integrations 32 .
Statistical and geospatial analyses. Collected data have been processed using: Exploratory Data Analysis (EDA) and Geostatistical Analysis (GA).
EDA was performed to evaluate the basic characteristics of the data and their statistical distribution by using numerical (i.e., calculation of summary statistics and statistical distribution of each variable) and graphical methods (i.e., histograms, box plots and normal probability plots). In particular, we used the Normal Probability Plots (NPP) to determine the occurrence of different overlapping geochemical populations and define threshold values by approximating linear segments on the point distribution in the graph 49  www.nature.com/scientificreports/ deviations or gaps in the NPP may indicate the presence of subpopulations (e.g., background and/or anomalies) separated by a threshold value [51][52][53] . GA was performed to understand and reconstruct the natural phenomena that govern the spatial behaviour of the studied variables. In particular, we used GA to visualize the samples distribution and the distribution of the measured values compared to their nearest neighbours, to study the spatial autocorrelation of the variables and elaborate a spatial model in order to estimate the variable values at unsampled locations and construct final prediction maps.
We accomplished this process according to the following steps: (1) Construction of experimental variograms of the studied variables. In particular experimental variograms and variogram surfaces were used to check the spatial continuity of the data distribution values and the presence of anisotropic phenomena (i.e., fault-related) acting along preferential directions; (2) Determination of the anisotropy (where present) which is important for defining parameters for the kriging estimation (i.e., directions and anisotropy ratio); (3) Construction of a spatial model, i.e., calculation of the main variogram parameters (nugget, range, sill) to be used in the kriging algorithm; (4) Preparation and validation of prediction maps (i.e., contour maps) by using ordinary kriging (OK) algorithm 25 .
Collected data were processed using the following software: ArcGIS Pro 2.7.0 (copyright 2020@Esri Inc.) and Surfer 23.1.162 (copyright 1993-2021, Golden Software, LLC) for the mapping process; Grapher 19.1.288 (copyright 1992-2021, Golden Software, LLC) and Statistica 12 (copyright Statsofts. Inc.) software for the numerical and graphical statistics.

RAI calculation.
To quantify the Tectonically Enhanced Radon (TER) we have applied the concept of Radon Activity Index (RAI). Over the years, many systems for classification radon concentration in soil gas have been elaborated to mathematically explain the geochemical activity of a fault zone. Seminsky et al. 11 proposed the relative index of radon activity K Q calculated as the ratio of the maximum 222 Rn concentration (Q max ) to the minimum 222 Rn concentration outside the fault zone (Q min ) for the classification of fault activity. Then, the values have been divided into five different levels of activity.
Inspired by this study, we have modified the concept of the index of radon activity, considering a Radon Activity Index (RAI) calculated as the ratio between the maximum 222 Rn value along the measured or estimated radon profiles to the background value (50 kBq·m −3 ) estimated over the area with the NPP method (see section geostatistical and spatial analysis in methods).
We have calculated the RAI considering only the radon peaks (radon values above the background threshold) relative to 21 (P4 to P24) estimated profiles, obtained by intersecting the radon grid map with 6 km long N-S sections (using Profile, Map Tools in Surfer 23.1.162) with a sampling step of 10 m and spaced 500 m from each other. Applying this method, each profile is composed of 600 estimated points of measurement ensuring good statistic representativeness. Based on the results obtained by RAI calculation, the values have been divided into four classes of geochemical activity: (i) RAI < 2: low activity, (ii) 2 < RAI < 3: medium activity; (iii) 3 < RAI < 6 high activity; (iv) RAI > 6 very high activity (see Supplementary Table S2).
The location of the RAI values is reported as classed post in a bivariate colour map (summary sketch map in Fig. 4) obtained using the bivariate colours option in ArcGIS Pro. The different colours in the bivariate map represent different values of the first (RAI) and the second (fault density) variable simultaneously.

Data availability
All data generated or analysed during this study are included in this article (and its Supplementary Information files).