Mapping radon-prone areas using γ-radiation dose rate and geological information

The national radon database currently has over 11 000 records of in-house radon measurements grouped by municipality. Several research groups have contributed to the database (Quindós et al 1991, Gutiérrez et al 1992, Amorós et al 1995, Baixeras et al 1996, Pérez et al 1996, Baeza et al 2003, Quindós et al 2011), the bulk of measurements having been collected by the University of Cantabria through different projects sponsored by CSN. The resulting sampling density is spatially heterogeneous, tending to be higher in regions with elevated radon level and at the most densely populated areas.

Roughly 90% of measurements were made with solid-state nuclear track detectors exposed for a period from three to six months, whereas the remaining 10% were obtained with charcoal canisters exposed over three days. Although, on a house-by-house basis, short-term measurements do not provide as reliable estimates of the annual radon level as the several-month-long measurements, an overall good agreement was found in surveys conducted in all municipalities where the two types of result could be compared (García-Talavera et al 2013). Furthermore, because our scope was to quantify risk rather than actual population exposure, we selected measurements made at ground or first floors.

The current radon database suffers from some limitations that hamper modelling approaches: it neither reports precise measurement location (only name of the village, town or city is provided) nor contains information on relevant dwelling characteristics (such as year of construction or type of foundation). Moreover, there are a limited number of measurements in comparison to the building stock. Out of the approximately 7900 municipalities in peninsular Spain, radon measurements have been made in 1170. Yet, for roughly one-third of them, the number of available measurements (n) is 1 or 2; 452 have 7 ≤ n < 26; and 25 have n ≥ 27.

Many countries have developed γ-radiation maps of their territories, as such maps are cost effective and have manifold applications: information on population radiation doses, documented baseline for accident and emergency situations, aid to raw materials or mineral exploration, etc.

The MARNA (Map of Natural Radiation of Spain) was undertaken as a collaborative project of CSN and ENUSA and joined later by several universities (Suárez et al 2000). Its main objective was to draw up a map of natural γ-radiation expressed as microroentgens per hour (μR h⁻¹) at 1 m above ground.

It combined information from various sources.

• Over 150 000 aerial measurements from uranium deposit exploration obtained in radiometric flights at an altitude of 120 m between 1967 and 1981. The detectors used were large volume NaI scintillators with an energy window from 0.4 to 2.8 MeV.
• Ground-based measurements, static and car-borne, also obtained in uranium exploration campaigns covering an energy window from 0.4 to 2.8 MeV.
• Over 450 000 ground-based static and car-borne measurements from project-specific campaigns using Saphimo SRAT and ES3 S equipment for total counting and an Exploranium GR-1300 for spectrometric measurements.

After correcting the airborne data for cosmic radiation, a good agreement was found between airborne and ground measurements (linear correlation coefficient greater than 0.92). To draw up the map to 1:1000 000 scale, mean radiometric values (corresponding to a nominal cell size of 7 km×5 km) were used. The 7 km×5 km grid of regularly spaced values was submitted to a kriging interpolation using Golden Software Surfer 7. In this step, an anisotropic variogram model was specified whose parameters were locally varied to describe the orientation of geological structures. The performance of the kriging model was verified by standard cross-validation.

Lithostratigraphic (LS) information was extracted from the 1:200 000-scale national map developed by the Institute of Geology and Mineral Exploration (IGME 2012a). This map represents cartographic units that are delimited mainly on the basis of lithic characteristics and stratigraphic position.

The metallogenic map at 1:200 000 scale (IGME 2012b) showing the nature and distribution of mineral resources was used to establish the location of uranium deposits and occurrences. The most important of the several tens of known uranium mineralised zones are located in the Schist–Greywacke Complex, a geological structure of the Central Iberian Hercynian belt. They mainly occur in fracture and breccia zones of Palaeozoic slates and granitic plutons.

If we denote by X_0.9 the 90th percentile of indoor radon concentration, X, in a given area, the problem of classifying it as radon prone or not can be formalised as a hypothesis test problem where the null hypothesis, H₀: X_0.9 ≥ 300 Bq m⁻³, is tested against the alternative hypothesis, H₁: X_0.9 < 300 Bq m⁻³. The critical point for rejecting H₀ is ${\tilde {X}}_{0.9}$, the upper bound on X_0.9, at the α-significance level.

As it is customary, we assume that indoor radon data, X, follow a lognormal distribution. We denote by m_g the sample geometric mean and by s_g the sample standard geometric deviation. Then, the random variable Y = logX is normally distributed with mean μ_y and standard deviation σ_y.

Given a random sample {x_i}, the upper bound on the 90th percentile of Y, Y_0.9, can be obtained as (Hahn and Meeker 2011)

$$\begin{eqnarray*} &&{\tilde {Y}}_{p}=\bar {y}+{g}_{(1-\alpha ;p;n)}^{\prime}{S}_{y}, \end{eqnarray*}$$

where $\bar {y}$ is the sample mean of {y_i} = {logx_i}, ${S}_{y}=\sqrt{n/(n-1)}{s}_{y}$, n is the sample size, s_y is the standard deviation of {y_i}, and the factors ${g}_{(\gamma ;p;n)}^{\prime}$ are given in table 1 of Odeh and Owen (1980).

Because of the relationship between the theoretical percentiles of X and Y, P{Y < Y_0.9} = 0.9↔P{X < e^Y_0.9} = 0.9, the upper bound on the 90th percentile of X is

$$\begin{eqnarray*} {\tilde {X}}_{0.9}&=&\exp \{ {\tilde {Y}}_{p}\} =\exp \{ \bar {y}+{g}_{(1-\alpha ;0.9 n)}^{\prime}{S}_{y}\} =\exp \{ \bar {y}\} \exp \{ {g}_{(1-\alpha ;0.9;n)}^{\prime}{S}_{y}\} \\ &=&{m}_{\mathrm{g}}{s}_{\mathrm{g}}\exp \left \{ {g}_{(1-\alpha ;0.9;n)}^{\prime}\sqrt{\frac{n}{n-1}}\right \} . \end{eqnarray*}$$

At the α-significance level, the null hypothesis can be rejected (thus concluding that a particular area is not radon prone) if and only if ${\tilde {X}}_{0.9}\lt 3 0 0~\mathrm{Bq}~{\mathrm{m}}^{-3}$. In particular, we set α = 0.10 (probability of type-I error), which allows for a 10% probability of failing to identify a true radon-prone area.

The above equation can be applied for n > 2, but small sample sizes would result in large probabilities of type II error (incorrect non-rejections of H₀). To control this effect, Faulkenberry and Weeks (1968), suggested that the sample size be large enough so that

• (i)  
  there is a large probability (1 − α) that the tolerance bound will be exceeded by at least 100p'% of the normal population, and
• (ii)  
  the probability (δ) that more than 100p*% of the population will exceed the said tolerance bound is small, where p' and p* are specified proportions, and p* > p'.

The probability of a successful demonstration is p_DEM = Pr(p > p') ≥ 1 − δ. For p' = 0.90; p* = 0.98 and δ = 0.10, we get a minimum sample size of 27. Consequently, we established this figure as the minimum sample size to apply the radon-prone hypothesis test.

Previous research (Cascón et al 2002) had shown that γ-dose rate across the country exhibits a linear relationship with ²²⁶Ra in soil, which, in turn, is known to be correlated with ²²²Rn indoors (Nason and Cohen 1980). This fact prompted a study on the usability of the γ-radiation map of Spain to predict radon levels directly in a northwestern granitic region (Quindós et al 2008). Three areas with different γ-radiation levels (<65–87 nGy h⁻¹, 65–130 nGy h⁻¹, and 87–>130 nGy h⁻¹) were measured for radon indoors. Based on the main statistical parameters for the three corresponding radon distributions, this study concluded that γ-radiation dose rate is a qualitative indicator of high indoor radon level rather than a good quantitative predictor.

We analysed the relationship between indoor radon concentration and γ-dose rate using data from the whole Spanish mainland. In particular, for each municipality, we considered the geometric mean (m_g) and the geometric standard deviation (s_g) of the measured radon concentrations, provided that there were more than six available measurements³. The corresponding average γ-dose rate was obtained from the MARNA map as follows: the built-up area of the municipality was determined using information provided by the National Centre for Geographic Information in a 1 km × 1 km grid. By a bilinear interpolation from the MARNA surfaces, the program Surfer 7 provided the γ-exposure value at the centre of every grid cell over which the built-up area spreads. We then computed the average of these values and converted it to absorbed dose in air using 1 R =0.008 76 Gy.

Our analysis showed that there is a statistically significant but rather weak correlation between radon geometric mean (GM) and average γ-dose rate (R_γ): Pearson r = 0.30, p = 0.001; Spearman ρ = 0.65, p = 0.000.

Radon geometric standard deviation (GSD), on the other hand, is not linearly correlated either to GM (r = 0.087; p = 0.082) or to R_γ (r =− 0.060; p = 0.229), although it exhibits a weak but significant rank correlation with the former variable (ρ = 0.108, p = 0.030). The distribution of GSD for all municipalities with n > 6 has a median value of 2.05, and a range going from 1.13 to 6.93.

The natural logarithm of GM is better correlated with R_γ (r = 0.61; p = 0.000). We performed a minimum least-squares analysis for R_γ versus ln(GM) (see figure 1), obtaining an R² value of 0.37. A very poor fit (R² = 0.05) was achieved when, instead of geometric means, 90th percentiles (X_0.9), assessed by the maximum likelihood parametric estimate, were considered.

Zoom In Zoom Out Reset image size

None of these regression models can explain a sufficient amount of variation to allow the construction of radon maps that carry an acceptable confidence of accuracy. Yet the relevant information to label an area as radon prone is not actually m_g nor X_0.9, but the probability of them exceeding a certain threshold value. This means that, to this purpose, radon indoors can be expressed as a categorical rather than as a quantitative variable, as presented in section 4.2.

First, we investigated the role of γ-dose rate as a categorical variable using the same dataset as in section 4.1: we selected cut-offs of 66 nGy h⁻¹ (∼7.5 μR h⁻¹) and 123 nGy h⁻¹ (∼14.0 μR h⁻¹) for the γ-dose rate because the categories defined by these values provide a good division of geological units (the >7.5 μR h⁻¹ surface roughly corresponds to the Pre-Mesozoic geological setting of peninsular Spain, whereas the 14.0 μR h⁻¹ contour line encloses some granitic outcrops and westernmost distinct metamorphic areas of the Iberian Massif—see figure 4).

We considered the distribution of radon geometric means for each group of municipalities defined according to the γ-dose rate categories established above (see box plot in figure 2). The three of them can be described by a lognormal distribution (Kolmogorov–Smirnov test: Z₀ = 0.514; p₀ = 0.954; Z₁ = 1.161; p₁ = 0.135; Z₂ = 1.089; p₂ = 0.186). Levene's test, performed on the log-transformed data, showed heterogeneous variances (W = 5.303; p = 0.005). Thus, we applied the Kruskal–Wallis test, which revealed that the three groups are statistically different (χ² = 151.4; p < 0.05).

Zoom In Zoom Out Reset image size

Next we carried out a correspondence analysis to explore the degree of association between the following categories of γ-dose rate and indoor radon (GM).

• <66 nGy h⁻¹ and <70 Bq m⁻³.
• 66–123 nGy h⁻¹ and 70–120 Bq m⁻³.
• >123 nGy h⁻¹ and >120 Bq m⁻³.

A good association was found for the three pairs (see figure 3(a)). Figure 3(b) shows a graphic representation of the associated contingency table. A geometric mean of 70 Bq m⁻³ corresponds (for n = 27 and s_g = 2.37, namely, the 75th percentile for the GSD distribution) to a 90% confidence upper bound on X_0.9 of approximately 300 Bq m⁻³. The 120 Bq m⁻³ cut-off was selected because it achieves the closest association to the >123 nGy h⁻¹ category.

Zoom In Zoom Out Reset image size

Out of 181 municipalities with average γ-dose rate < 66 nGy h⁻¹, nearly 90% (95% CI=(0.85,0.94) estimated by the modified Wald method; Agresti and Coull 1998) have radon geometric means lower than 70 Bq m⁻³. A higher proportion, 97% (95% CI=(0.93,0.99), likewise calculated), have 90th percentiles lower than 300 Bq m⁻³.

Therefore, selecting a cut-off point of 66 nGy h⁻¹ for a screening test based on γ-dose rate can be a valid option to discriminate non-radon-prone municipalities. We considered two such screening tests using the following criteria, respectively, to define radon-prone areas: (i) m_g ≥ 70 Bq m⁻³; (ii) X_0.9 ≥ 300 Bq m⁻³. Table 1 shows their positive predictive value (PPV, probability that a municipality is radon prone given that the test result is positive) and their negative predictive value (NPV, probability that the municipality is non-radon-prone given that the test result is negative). Both tests have a poor PPV but achieve a high NPV, notably for X_0.9.

Because of the high NPV, regardless of the selected point estimate, all areas with γ-dose rate lower than 66 nGy h⁻¹ (∼7.5 μR h⁻¹) were directly classified as non-radon-prone—see map in figure 4. Most of these correspond to Tertiary and Quaternary soils developed over calcareous sedimentary rocks.

Zoom In Zoom Out Reset image size

Conversely, for areas with γ-dose rate in excess of 66 nGy h⁻¹, γ-radiation does not exhibit a strong enough control on radon indoors to make a reliable classification (see the first two bars in figure 3(b)). Even the highest γ-dose category contains a significant proportion (∼20%) of non-radon-prone municipalities.

The control of geology on indoor radon levels is well recognised (Miles and Ball 1996, Mikšová and Barnet 2002, Kemski et al 2002, Appleton and Miles 2010, Friedmann and Gröller 2010). Consequently, for those parts of the country with γ-dose rate in excess of 66 nGy h⁻¹ (i.e., >7.5 μR h⁻¹ surface from map in figure 4), we included geological information in the analysis in order to establish radon-prone areas at geological-formation level rather than at municipality level. Specifically, we considered the lithostratigraphic (LS) map of Spain.

To analyse radon data, we selected all LS units that intersected in more than 10% of their surface with populated regions exhibiting γ-dose rate >66 nGy h⁻¹. Some of these units were further regrouped based on composition, age and permeability (see table 2). The LS unit has, on the one hand, a characteristic range of ²²⁶Ra concentration and, on the other hand, a typical permeability⁴ value, both of which have a link to indoor radon. Yet, even within the same LS type, a substantial variability of indoor radon levels could be found due to the influence of unaccounted-for variables (such as climate, dwelling characteristics, or depth to the water table).

Table 2.  Statistical parameters of radon levels and average γ-dose rate (${\hat {R}}_{\gamma }$) in lithostratigraphic units within regions of γ-dose rate values predominantly higher than 66 nGy h⁻¹. Those whose upper bound on the 90th percentile (${\tilde {X}}_{0.9}$) of measured radon concentrations exceeds 300 Bq m⁻³ are classified as radon prone.

Unit description	Sample size	${\hat {R}}_{\gamma }$ (nGy h−1)	mg(Bq m−3)	sg	${\tilde {X}}_{0.9}$ (Bq m−3)
Hercynian acid plutonic rocks (Palaeozoic)	770	131	118.0	2.7	452.5
Hercynian arkosic sandstone and schist (Palaeozoic)	96	121	57.6	2.8	261.9
Sedimentary deposits of Tertiary age (Guadiana basin)	164	117	55.7	1.8	130.9
Slates, metagraywacke and mica-schist (Precambrian–Palaeozoic)	331	114	98.0	2.7	375.4
Metamorphic rocks (Precambrian–Cambrian–Ordovician)	85	110	131.9	2.3	459.6
Tertiary siliciclastic deposits, mainly arkosic sandstone and mudstone (Tagus–intra-Betic–Guadalquivir)	131	110	70.1	2.0	185.8
Precambrian slates and greywacke (Schist–Greywacke Complex)
SC	119	107	105.1	2.9	493.9
SCb	94		53.1	2.0	150.6
Conglomerate/lutite/coal (Palaeozoic)	31	96	83.1	1.7	203.2
Olot volcanic rocks (Quaternary)	68	83	62.1	2.3	218.1
Palaeozoic microconglomerate, sandstone and shale (Cantabrian and Western Astur–Leonian Zones)	82	80	70.7	2.0	207.7
Alluvial basin-fill deposits of Tertiary age (Douro basin)	95	78	44.5	2.1	130.5
Quaternary basin-fill deposits along stream flood plains, salt marshes, estuaries, etc	491	75	43.6	2.5	152.7
Flysch sandstone and lutite alternations (Ebro basin and Pyrenees)	45	68	44.7	2.3	170.6
Alluvial basin-fill deposits of Tertiary age (Guadalquivir basin)	60	64	30.0	2.6	129.5

We first analysed whether radon variability within every LS unit is a random fluctuation about the mean or conforms to a spatially structured pattern. In the former case, the radon-prone hypothesis test could be directly applied on the combined raw data. The latter case would possibly require a different approach, since spatial patterns in the data could impair the performance of the standard statistical test of a hypothesis (Legendre 1993).

Because no geo-referenced information on individual radon measurements was available, we examined the variability of aggregate data expressed as the municipalities' geometric means (considering only those calculated from more than six measurements). Besides, to avoid bias in the results due to lack of precise measurement location, we excluded municipalities whose built-up area spreads over more than one LS unit. Within every LS unit, data are consistent with a lognormal distribution, as verified by the Shapiro–Wilk test at p = 0.05. We thus obtained the log-transformed geometric means in order to achieve normally distributed attribute data.

Unfortunately, in most LS units, the small number of municipalities, together with an irregular spatial configuration, made it impossible to conduct reliable spatial statistical analyses (it is, for instance, difficult to properly estimate a variogram with fewer than 50 data points; Olea 2006). Visual inspection did not reveal any global pattern except for slates and schist (SC) of the Schist–Greywacke Complex (see figure 5). As for the other LS units, the lack of spatial structure may be statistically confirmed on the most extensive ones (namely Hercynian acid plutonic (AP) rocks and Quaternary (Q) sedimentary material) for which sample sizes are large enough (>50 municipalities).

Zoom In Zoom Out Reset image size

We obtained the experimental variograms for units AP and Q (figure 6). Both of them correspond to a pure nugget effect, as can be seen from the theoretical model fit. This model entails a complete lack of auto-correlation at distances larger than the lag bin size.

Zoom In Zoom Out Reset image size

The same results were confirmed by means of Moran's I (Cliff and Ord 1973). Moran's I ranges between −1 (perfect negative auto-correlation at global level) and 1 (perfect positive global auto-correlation), with values close to 0 indicating a random distribution. We used an inverse-distance weight matrix to compute the statistic, and a Monte Carlo permutation test (i = 1000) to assess its statistical significance. For plutonic rocks, we obtained I_obs = 0.066 (p = 0.051), whereas for Quaternary materials I_obs = 0.039 (p = 0.142). Both I_obs and p values indicate a lack of global auto-correlation.

The lack of a global spatial pattern (gradients or large-scale patches) suggests that precipitation and temperature regimes, which exhibit marked gradients N–S/SE through these two extensive LS units, do not exert a noticeable influence on radon indoors. Furthermore, given that the two units exhibit opposite behaviour regarding radon transport (impermeable matrix with highly permeable fracture zones versus highly permeable homogeneous material), it is reasonable to assume that this result may be extrapolated to all LS units excepting unit SC, which needs further analysis.

There are two issues regarding radon data from unit SC. First, one of the sampled municipalities displays a clearly anomalous behaviour (see box plot in figure 5) and was thus excluded from the study. The very high radon values (up to 15 000 Bq m⁻³) measured in this municipality can be attributed to the existence of an underlying uranium deposit (Arribas et al 1983).

Although this is a unique situation, other municipalities have close uranium deposits that could potentially increase indoor radon levels. The diffusion length of this gas in soil is only about 1.2 m, but advective transport by carrier gas also occurs at long distances (several hundred metres) (Etiope and Martinelli 2002). In addition, surface soils in the proximity of uranium deposits show enhanced ²²⁶Ra concentrations (up to two orders of magnitude higher than background levels), although this is a very local effect (Baeza et al 1994). For unit SC, we compared indoor radon data {x_i} from five municipalities near uranium deposits (0.3–2 km) with those from the rest of the locations. The F test indicates equal variances for both groups (F = 0.5821; p = 0.094), whereas the two sample t-test proves equal means (t = 0.5346; p = 0.5929). Similar results are actually obtained when the analysis is replicated for unit AP, which also hosts several uranium deposits: F = 0.5841 (p = 0.062); t = 0.6405 (p = 0.522). Based on these results, we concluded that uranium deposits in granite and schistose rock do not have an influence on indoor radon levels at nearby population centres, at least when located at distances larger than 0.3 km (from the built-up area to the boundary of the ore deposit).

The second issue regarding radon levels in unit SC is the presence of a visually apparent large patch of minima that can be observed in figure 5. Most of the municipalities forming this patch are located in an isolated mountain range that forms a series of alternating ridgelines and valleys. They were considered as a distinct LS unit (SCb) solely in order to detrend the data and prevent bias in statistical inference.

We then grouped the individual ${\mathop{\sum }\nolimits }_{j=1}^{{k}_{m}}{\mathop{\sum }\nolimits }_{i=1}^{{n}_{j}}{n}_{i,j}$ radon measurements made in dwellings at the k_m municipalities in each LS unit, m. For each of these data samples, {x_i}_m, we performed a normality test on the log-transformed values. If the sample size was larger than 50, we used the Kolmogorov–Smirnov test. Otherwise, we applied the Shapiro–Wilk test. All data sets follow lognormal distributions and, thereby, the statistical hypothesis test described in section 3.1 was performed. Table 2 shows ${\tilde {X}}_{0.9}$ for the different LS units, as well as other relevant statistical parameters.

The resulting map of currently identified radon-prone areas is shown in figure 7. All units classified as radon prone have average γ-dose rates greater than 100 nGy h⁻¹, whereas some units with average values as high as 120 nGy h⁻¹ are not radon prone. Also, at the geological unit level, γ-dose rate plays a role in determining radon risk but is not a good indicator above a certain threshold.

Zoom In Zoom Out Reset image size