Comparison of radon mapping methods for the delineation of radon priority areas – an exercise

Background: Many different methods are applied for radon mapping depending on the purpose of the map and the data that are available. In addition, the definitions of radon priority areas (RPA) in EU Member States, as requested in the new European EURATOM BSS (1), are diverse. Objective: 1) Comparison of methods for mapping geogenic and indoor radon, 2) the possible transferability of a mapping method developed in one region to other regions and 3) the evaluation of the impact of different mapping methods on the delineation of RPAs. Design: Different mapping methods and several RPA definitions were applied to the same data sets from six municipalities in Austria and Cantabria, Spain. Results: Some mapping methods revealed a satisfying degree of agreement, but relevant differences were also observed. The chosen threshold for RPA classification has a major impact, depending on the level of radon concentration in the area. The resulting maps were compared regarding the spatial estimates and the delinea‐ tion of RPAs. Conclusions: Not every mapping method is suitable for every available data set. Data robustness and harmon‐ isation are the main requirements, especially if the used data set is not designed for a specific technique. Different mapping methods often deliver similar results in RPA classification. The definition of thresholds for the classification and delineation of RPAs is a guidance factor in the mapping process and is as relevant as harmonising mapping methods depending on the radon levels in the area.

The definition of RPA in the EU-BSS allows a wide range of interpretation, and therefore, different concepts and methodologies have been proposed and already adopted in some countries (2,3). Radon maps have existed in several countries for many years as part of national radon strategies even before the new EU-BSS became effective (4)(5)(6)(7)(8)(9)(10). The applied mapping methods and the visualisation are very different among the various countries, depending on the purpose of the map and the data available (11). These methods are based on different developments, strategies and ideas in radon protection for many years in the countries, and most of the time, the basic mapping strategies and methods applied in a country remain unchanged, even when revised or new legal requirements were applied. Consequently, a basic bottom-up harmonisation approach (same methodology everywhere) of mapping methods or definition of RPA will not be enforceable. Therefore, comparison, evaluation and discussion for possibilities of top-down harmonisation (different methodologies are normalised to common standards) are important. Harmonised evaluation, classification and display of the radon potential are important for better comparability and compatibility between regions or countries and should serve as a basis for appropriate and consistent radon protection measures for the population. Within the European research project 'Metrology for Radon Monitoring (MetroRADON)' (12), the goal was to develop reliable techniques and methodologies to enable SI traceable radon measurements and calibrations at low radon concentrations, including also the task of harmonisation of radon data and RPA. The aim was to develop a strategy to harmonise defined RPA across borders. In this framework, studies and exercises were carried out based on literature, available data and case studies as a basis to evaluate the situation and develop strategies. Results of these studies are reported in the MetroRADON deliverable (13) and will also be discussed in journal articles, for example, the topic of causes and effects of lack of compatibility between maps and possible methods for harmonisation (14). In this article, the results of the 'radon mapping exercise' carried out within MetroRADON are presented and discussed. The idea of the exercise was to evaluate if available and established mapping methods can be applied to a data set from another area and if different mapping methods applied on the same data deliver comparable results. So, different radon mapping methods already used in countries for RPA definitions were applied to two harmonised data sets of various variables (e.g. indoor radon, gamma dose rate, geology, soil gas radon [SGR]) to evaluate the impact of the different techniques on the delineation of RPAs, as well as their potential applications to other countries.

Data set description
The basis for the realisation of the mapping exercise was the availability of suitable data sets. The demands were: to have more than one data set, preferably from different countries or regions; the data set should include various variables which could be interesting for mapping (e.g. information about indoor radon, SGR, geology, geogenic parameters); and to have the permission that the data set can be shared with the participants of the exercise and the results used for the project. Based on these requirements, the selected data sets for the mapping exercise are from different radon measurement campaigns in six municipalities in Austria and in Cantabria, Spain.
The data include indoor radon concentration (IRC) measurements in dwellings, building characteristics of measured dwellings, SGR concentration, soil permeability, radionuclide concentrations ( 226 Ra, 228 Ra, 210 Pb, 228 Th,   232 Th, 238 U, 40 K) in soil samples, ambient dose rate (ADR) and maps of geology, soil type and airborne radiometry (see Table 1). All data are georeferenced and provided in shape files or TIFF raster files.
The Austrian data set covers six municipalities and is separated into two distinct areas located in the North and in the South of Austria (AUT North and AUT South, respectively, Fig. 1, bottom). Each area consists of three adjacent municipalities with an overall area of about 220 km² (40 km² in AUT North, 180 km² in AUT South).  Most of the data were collected during detailed measurement campaigns carried out between 2010 and 2012. The survey data are supplemented with data obtained from literature (see Table 1). The area AUT North is located in the Bohemian Massif which is characterised by high geogenic radon potential (GRP) due to the predominant presence of granites and gneiss outcrops. It shows homogeneous geological features with a granitic pluton and interlaying migmatites (metamorphic rocks with granitic parent rock). The geology of AUT South is comparatively heterogeneous and is characterised by a variety of felsic igneous and metamorphic rocks with high radon potential, but also sedimentary units with low radon potential.
References to the geological maps are given in Table 1.
More details about the radon survey and data are discussed in Refs. (15)(16)(17)(18)(19) and more information about radon and geology in Austria are given in Ref. (10,20). The Spanish data set covers the region of Cantabria ( Fig. 1, top) having a total area of about 5,300 km². The data set consists of different measurements of IRC, SGR, ADR and data compiled from literature (see Table 1). The geology of Cantabria is mainly characterised by detritic sediments and carbonate rocks, which usually show low to intermediate radon potential; however, the high permeability of the fractured carbonates can result in locally higher GRP. The metasediments located in the western part of the region and local volcanoclastic formations usually show a low GRP. In general, compared to the Austrian regions, Cantabria has lower GRP. References to the geological maps are given in Table 1. More details about radon mapping in Spain are discussed in (21)(22)(23)(24). Figure 1 shows the studied areas in Austria and Spain. The figure also includes classed post maps of the the IRC measurements and examples of other available data layers (e.g. geology, water conditions and lithology).
The data sets differ in basic characteristics such as size, sample density, quality and resolution. The characterisation of RPA for the two data sets may require adequate data preparation according to the different mapping methods. Table 1 reports an overview and a comparison of the Austrian and the Cantabrian data sets regarding data density, similarity, source (e.g. measured or derived from literature) and number of measurements (where applicable).
Detailed analyses were carried out for all variables of all data sets (AUT North, AUT South, Cantabria). Descriptive statistical analysis was performed by box-plot graphs and by other statistical methods (Kruskal-Wallis test, Spearman's rank correlation, variograms). Details can be found in the final report of the exercise (35). Table 2 shows a qualitative summary of the statistical and spatial correlations of the analysed variables. The analysis of the datasets indicates that the different regions do not show the same results regarding the correlation between quantities, differences of a quantity within groups (e.g. bedrock type, soil grain size, water content in the soil, permeability, building characteristics) and strength of spatial correlations.

Methods
Different mapping methods are discussed within the exercise. The idea was to include as many mapping methods as possible in the exercise, which were already used for radon mapping in countries or were suggested by experts. Therefore, experts from different countries were invited to participate in the exercise and apply their respective mapping method to the provided exercise data sets. Not all invited experts had the time or resources to perform the exercise (as we could not provide funding for the external participants for this work within the project). The mapping methods and results of those who agreed to participate in the exercise were included in this article. In addition, basic statistics of indoor radon data was performed, as basic statistic methods are also used for radon mapping in several countries.
In the following, these methods are briefly described, and examples of comparisons are examined in the Results section. More details about the methods and results can be found in the final report of the exercise within the MetroRADON project (35) and the reported specific literature.

Basic Statistics of Indoor Radon Data
The definition of RPA by using IRC data commonly follows one of two basic concepts: 1) the mean IRC (e.g. AM, GM) of the area is compared to a threshold (e.g. 300 Bq/m³) and 2) the percentage of measurements exceeding a threshold (RL) in an area is compared to a percentage threshold (e.g. 10%). Common approaches to define RPA use IRC thresholds (RLs) ranging from 100 to 300 Bq/m³ and percentage thresholds ranging from 1% to 30%. Descriptive statistical analysis was performed for the exercise data in the light of these two basic RPA concepts, as a common mapping method. Results are shown in Table 3 and Fig. 2.
Generalised Additive Mixed Model (GAMM) GAMM (36) is used to estimate the IRC (as a dependent variable) by using the correlation with some explanatory variables. The method is based on Ref. (37). For Austria, the additive mixed model: 1,..., 1,..., ) ( 0, ) is fitted to the data set, whereby the living unit u j is taken as a random effect, thus introducing a positive correlation Comparison of radon mapping methods  of measurements within the same living unit, because u j is common to all measurements within unit (j). A slightly different additive model without random effects is used for the Cantabrian data set, because measurements from a relatively large area are assigned to a particular location.
Influencing factors, such as geology, in such an area could be inherently different, which would contradict the positive correlation induced by the random effect. In both cases, the smooth functions s(.) pertain to the class of thin plate regression splines. The z ij terms represent explanatory variables and the pair (x j , y j ) represents the coordinates of a living unit or location j. The final model should only contain variables that show a significant influence on log(IRC).
To identify these variables, a stepwise forward selection using a 5-fold random cross validation was applied.
Variables with the highest explanatory power were chosen for the final model. Non-relevant variables result in non-significant improvements in cross validation error. For Cantabria, the following explanatory variables based on the fivefold random cross validation were used for the final model: soil-type, ADR, K 2 O and Th content in soil. For the Austrian data sets, the following building characteristics were used: room earthbound (yes/no), floor, type of walls, type of basement (yes/no/partly), solitary building (yes/no) and type of bedrock. Additionally, the type of foundation and the U content were selected for the AUT North data set, while the fraction of measurement time in winter, the number of dwelling units, tightness of windows and water saturation of the soil were selected for the AUT South data set. The final models are fitted using these variables to predict log(IRC) according to specified grids (10 × 10 km for Cantabria and 2 × 2 km for Austria). Some results can be found in Table 3.

Ordinary Kriging (OK) and Indicator Kriging (IK)
The kriging method (38)(39)(40) was performed for the prediction of IRC in both Austria and Cantabria. For Austria, an analysis of variance (ANOVA) was performed for the following target variables: IRC (only ground floor measurements were considered), SGR, soil permeability and the GRP after (41) (as function of SGR and soil permeability). ANOVA revealed significant (P < 0.05) differences for the target variables dependent on pedology and geology. Considering the high density of IRC measurements in populated areas, a pure geostatistical approach using OK and IK without any additional predictor seemed to be sufficient to estimate the radon risk for populated areas. First, the spatial autocorrelation of IRC was tested by calculating variograms. Based on these variogram models and the empirical data, IRC was kriged for a raster cell size of 200 m. Due to the low range of spatial autocorrelation, the estimates at large distances from the nearest observation (> 1 km) are equivalent to the mean of the whole area. The radon risk mapping was conducted using IK. For this purpose, IRC was transformed into a binary code with 0 for all observations that are smaller than 300 Bq/m³ and 1 for all observations that are greater or equal to 300 Bq/m³. Another variogram model was fitted to the binary coded data.
As the IRC data from Cantabria have no exact coordinates and no information was available about the floor of the building in which the measurement was performed, a different method was used. To make the data ready for kriging, all measurements from one municipality were merged into one value by calculating the arithmetic mean (AM). Thus, each unique location is assigned to one value for IRC. Spatial autocorrelation of IRC was tested but not detected, that is, the empirical data could not be fitted in a meaningful way to the variogram model. Hence, kriging of IRC was not a feasible option for the delineation of radon risk areas in Cantabria. Instead of a geostatistical analysis of IRC, the GRP was calculated as a function of SGR and soil permeability. Kriging, based on an exponential variogram model, was conducted for SGR.

Comparison of radon mapping methods
Data on soil permeability was assigned to five permeability classes, depending on the lithological type. Consequently, it was not possible to use the Neznal GRP (41). According to the Cantabrian data set, the GRP was defined as: Examples of the results by Kriging methods are shown in Fig. 4 (left hand side) and Fig. 7 (right hand side).

Empirical Bayesian Kriging Regression (EBKR)
EBKR is a geostatistical interpolation method that combines kriging interpolation and ordinary least square regression providing an accurate prediction of moderately non-stationary data at a local scale. It uses a dependent variable measured at point locations and known potentially correlated explanatory variables, as raster grids (42,43). EBKR estimates multiple semivariogram models instead of a single variogram by repeated simulations, thus accounting for the uncertainty introduced in the calculation of variogram parameters and providing a better accuracy than other kriging techniques. In common kriging methods, the prediction at unknown locations considers the nearby known data. This may result in the underestimation of the prediction standard errors caused by the uncertainty of semivariogram parameters. In contrast, EBK uses an intrinsic random function as the kriging model, differently than the other kriging methods, and does not assume a tendency toward the overall mean; this results in the same probability for large deviations to get larger or smaller. Furthermore, EBKR also considers the presence of a multicollinearity among the explanatory variables by using the principal components in the regression model. Each principal component captures a certain proportion of the total variability (set to 75%) of the explanatory variables. Cross-validation method was used to estimate the performance of the model.
In this work, the estimation by EBKR uses radon concentration in soil gas as response variable and raster layers of permeability, ADR, K-40, U-238, Th-232, fault density, presence of karst areas as predictors with a resolution of 500 × 500 m. Figure 4 (right hand side) shows the result of EBKR for GRP mapping of Cantabria.

Belgian Radon Risk Mapping Software (BRRMS)
Cinelli et al. (44) developed the method and the corresponding software is described in Ref. (45). The principle is to map the variations of the radon risk within geological units with the moving average method, while geological units with significantly different levels of risk are considered separately. When contiguous geological units have similar mean IRC levels, they are treated as a single unit. Within a given unit, the moving average of the nearest 20 data points is calculated (more precisely, the log mean, or the log median) for any chosen coordinate set, for example, the nodes of a square grid. The percentage of data locally exceeding a chosen threshold is also predicted, assuming lognormal distribution. The threshold used in the exercise is the European reference level of 300 Bq/m³, and the lognormal distribution is only fitted to data above the median (46). The method does not include a classification of the nodes. A classification in five risk classes is used in the Belgian method for municipalities (47), but was not included in the software.
Only the Austrian data set was used for this method, and only the highest concentration, measured on the ground floor, was kept for each living unit. The Austrian data sets are from two distinct rather small radon-affected areas. Each area includes different geological formations. However, the radon statistics give rather similar values for the geometrical mean (GM) IRC in the different geological units of each area, and therefore, they were considered as a single mapping unit. In Fig. 7 (left hand side), one example of the results is shown.

Results and discussion
The data sets for the exercise are complex, and correlations between variables were less significant than expected ( Table 2). The Austrian data sets represent only small areas (6 municipalities), which seems to be too small and geologically homogenous with respect to radon risk for geogenic correlations and modelling. The Cantabrian data set represents a larger area, but the data came from different surveys and literature (see Table 1). In addition, the Cantabrian data set has low sampling density and no exact coordinates for IRC, which makes the use of IRC for modelling challenging.
The fact that the data are inhomogeneous and not perfect in several aspects makes the exercise realistic, since, in practice, most of the time the available data for mapping are neither as perfect nor as complete as would be desirable. Consequently, the exercise shows how mapping methods can perform with incomplete or heterogeneous data sets, and how the classification of RPA can be done with them.
The data sets required adequate data manipulations to apply the different mapping methods, and not all data were used for each mapping method. Further, certain data characteristics (e.g. different input variables, not all variables which are needed for a method exist in all data sets) inhibit the application of each mapping method for every dataset. Table 4 gives an overview of the applied mapping methods and the data that were used for the respective method. In general, mapping methods are mostly specified to use either IRC or geogenic variables as input variables. The BRRMS combines IRC and geogenic variables by considering geological units. The methods using IRC with building characteristics could only be applied for the Austrian data sets, as no information about building characteristics is included in the Cantabrian data set. Only the GAMM method used all available variables for the Austria and the Cantabrian data sets and selected the relevant explanatory variables in a stepwise forward method. Apart from the basic statistic methods, all applied methods used interpolations to map the radon concentration, radon potential or the radon risk.
A summary of the results for Cantabria and the six municipalities in Austria for IRC derived from the different methods is shown in Table 3. The table gives an arithmetic/geometric mean/median value for the IRC or the percentage of measurements above 300 Bq/m³ in Cantabria and each of the six Austrian municipalities (Mun.). The methods which delivered results for grid cells were aggregated for the region of Cantabria and the municipalities in Austria. For this purpose, all grid cells were used, and the target variable of the corresponding method (median, geometric mean, arithmetic mean) was used as aggregate. This simple approach was chosen only to give an overview of the results derived from different methods based on administrative areas (province, municipalities), which RPA delineation is mostly based on. The results show that the predicted radon concentration is clearly lower for all methods in Cantabria than in Austria (see also Fig. 2) and also lower in the three municipalities in AUT South compared to AUT North. The GM of Cantabria data from basic statistics and the GAMM correspond very well, also for AUT Mun. 2 and 4. For the other municipalities, it deviates quite strongly, especially for Mun. 5 and 6. The BRRMS median (Med) concentration per municipality compared to basic statistics median deviates about 10-30%, and the deviation is stronger for the values of percentages above 300 Bq/m³. The OK IRC and BRRMS IRC predictions per municipality deliver higher values than the basic statistics, except for Mun. 3.
In the following, we will describe the consistency/comparability of different methods through a few examples.
Correlation analysis was performed only for methods, which provided the same variable as result (IRC, GRP) and the results were aggregated to the same grid. Figure 3 compares the GRP predictions for Cantabria obtained by applying EBKR and OK. The data were aggregated into a 5 × 5 km grid and the coefficient of determination (R 2 ) is 0.59. The correlation between the two methods for the area is good (acceptable). The GRP predictions of the two methods are displayed in the map in Fig. 4. The two maps show a corresponding pattern, with only some higher GRP in the North of Cantabria. OK predictions are generally higher for the central and southern part of Cantabria (Fig. 5). Possible reasons for this observation are seen as 1) utilisation of co-variable data by EBKR but not by OK and 2) predictions by OK beyond the range of spatial auto-correlation. In more  Comparison of radon mapping methods detail, for 1) OK uses only nearby observations for predictions whereas EBKR considers many environmental co-variables (such as geology). If, for instance, a geological unit with a small spatial extent and medium/too high risk exists within an area with general low risk, for OK a single medium/high Rn measurement would affect the whole area within the range of spatial auto-correlation by increasing the estimate irrespective of the environmental  setting. In contrast, for EBKR higher predictions would be more restricted to the respective geological unit where the medium/high Rn concentration was observed. For 2) low sampling density in certain areas could result in some cells being located beyond the range of spatial auto-correlation. This would be especially problematic for OK, because for these cells, the OK estimate tends towards the mean of the observational data. The mean of the observational data might in turn be influenced by more measurements from high Rn concentration areas (in case of Cantabria more measurements were conducted in the north where Rn concentration is higher). A detailed evaluation of performance of the individual maps in terms of its actual accuracy would have required independent test data that are currently not available or exhaustive simulation (e.g. via repeated cross-validation), which was beyond the scope of this study. However, performance assessment for different mapping methods applied for radon risk evaluation is certainly an important task that deserves a more detailed analysis in the future. Figure 6 compares the BRRMS method with the IK method for the predicted percentage of measurements above 300 Bq/m³ for the area AUT North. As basis for comparison, the coarser 500 × 500 m grid of the BRRMS was used and compared with the cell of the 200 × 200 m kriging raster closest to the midpoint of the BRRMS grid cell. The coefficient of determination (R 2 ) is 0.41, which is still a satisfying correlation. In Figure 7, the results of the two methods are displayed as maps. The two maps are similar, showing the highest radon potential in the centre. In general, the predicted IRC by the BRRMS method is higher than the one by IK.
Finally, we evaluated how the different results provided by different mapping methods would have an impact on the classification or delineation of RPAs. As discussed above, different definitions of RPA concepts were chosen   for the different methods for the same municipality impact the RPA classification, when the threshold is chosen in the range of the variability of the results (e.g. 300 Bq/m³, as shown in the example). If the threshold is set to 100 Bq/m³, all municipalities are classified the same, as this threshold does not lie within the range of the measurement/prediction results, and therefore, the variability of the results among the different methods does not impact the classification of RPAs.
If the threshold of percentage of measurements/predictions is set to 30% (over 300 Bq/m³), all municipalities in AUT North would be classified as RPAs with all three applied methods, as well as all six municipalities for the BRRMS method. Applying the commonly used definition of RPA in Europe (10% of dwellings above 300 Bq/m³), all six municipalities in Austria would clearly be considered as RPAs, independent of the mapping method. As discussed above, the variability of the results of the different methods only impacts the classification of RPA when the set threshold lies within the range of the predicted/measured results. Figure 8 displays the results for the three municipalities of the AUT North and Austria South areas for the different methods. The results (AM/GM/Med) per municipality for the respective methods are plotted, and the colouring shows for which threshold the municipality would be considered to be RPA (orange) and non-RPA (grey).

Conclusions
The evaluation of IRC and GRP in Europe and their harmonisation between countries and across borders constituted one of the main objectives of the MetroRADON project. This exercise-work, conducted within a specific task of the project, was aimed at the definitions of RPAs by using different mapping techniques. Results highlight that the application of a mapping method using data sets, not designed for the specific requirements of the used mapping method, is challenging. Usually, data sets always have specific characteristics and are rarely comparable, even for the same variable. Therefore, harmonisation is a challenge. However, some of the mapping methods used in this exercise show quite good correlations for the predicted cells, indicating that the used methods should in principle be interchangeable for harmonisation purposes. In general, the selection of a mapping method for a certain area will strongly depend on the available data sets and their statistical properties. Therefore, not all mapping methods are usable for all data sets or areas, depending especially on data quality, sampling density or natural heterogeneity of the mapping area. For harmonisation purposes of mapping at a European scale, a method using less parameters might be preferable, as it would be easier to apply to different data sets.
If a survey for delineation of RPA is started from scratch in a country, the mapping and display/classification methods (e.g. % above RL in administrative area) should be decided at the beginning, so that the survey design (e.g. sampling density and analysed parameters) can be optimised to these requirements.
Usually, the final goal of radon mapping is the delineation of RPAs, as this is requested in the EU-BSS. It was shown in this exercise that independent of the applied method for large intervals of classification thresholds, the same RPA classification is predicted. Different methods often deliver same results in RPA classification, according to the definition of RPAs. Problems emerge if the classification thresholds are close to mean IRC levels. In this case, small differences in estimated IRC between mapping methods can impact the RPA identification of the study area. The definition of the threshold values is an essential factor in the process of delineation of RPA. The definition of RPA is in general the most important factor that contributes to disharmony between RPA maps, and its harmonisation is as relevant as harmonising mapping methods.