Bayesian spatial modelling of terrestrial radiation in Switzerland

The geographic variation of terrestrial radiation can be exploited in epidemiological studies of the health effects of protracted low-dose exposure. Various methods have been applied to derive maps of this variation. We aimed to construct a map of terrestrial radiation for Switzerland. We used airborne $\gamma$-spectrometry measurements to model the ambient dose rates from terrestrial radiation through a Bayesian mixed-effects model and conducted inference using Integrated Nested Laplace Approximation (INLA). We predicted higher levels of ambient dose rates in the alpine regions and Ticino compared with the western and northern parts of Switzerland. We provide a map that can be used for exposure assessment in epidemiological studies and as a baseline map for assessing potential contamination.


Introduction
Terrestrial radiation stems from radio nuclei contained in the topsoil. The main contribution comes from 40 K and the elements of the Uranium and Thorium decay chains. The levels of ambient dose rates from terrestrial origin vary in space, depending on the local concentrations of the respective radioactive isotopes in the soil (UNSCEAR, 2018).
This geographic variation has been exploited in epidemiological studies on the health effects of protracted low dose exposures. Recent work has specifically looked at possible links between risk of childhood cancer and background ionising radiation (Mazzei-Abba et al., 2020;Demoury et al., 2017;Spix et al., 2017;Kendall et al., 2013). The effects on cancer risks of protracted low-dose ionising radiation are expected to be small and require large sample sizes to be detected. As direct measurements of doses are not feasible for large numbers of study participants, these studies assessed exposures using geographic exposure models and geocoded residential address information.
A variety of methods have been used to develop geographical models for terrestrial background radiation for exposure assessment. In Great Britain, an ordinary least squares (OLS) regression including predictors such as bedrock classes and radiation measurement means by county provided the best performance (Chernyavskiy et al., 2016;Kendall et al., 2016) among various modelling approaches. In France, Warnery et al. (2015) applied ordinary kriging and multi-collocated co-kriging to a large dataset of dosimetry measurements conducted in veterinary clinics. The two methods performed similarly. However, the co-kriging approach, which jointly modelled geogenic uranium potential with terrestrial radiation, showed more detailed spatial features. In Switzerland, Rybach et al. (1997, 2002 used inverse distance interpolation to derive maps of terrestrial gamma radiation from naturally occurring radionuclides and from 137 Cs (fallout from the Chernobyl nuclear accident) based on a heterogeneous set of measurements including in situ measurements, laboratory measurements of stone samples and airborne γ-spectrometry. These estimates were combined with dose rates from cosmic radiation, calculated as a function of elevation, to obtain a map of total external background radiation.
Integrated nested Laplace approximations (INLA) used together with stochastic partial differential equations (SPDE) allow fitting models involving Gaussian random fields (GRF) to large data sets at reasonable computation costs by establishing an explicit link between the GRF and Gaussian Markov random fields (GMRF). This method has for instance been applied to the spatial prediction of soil pH and elemental concentrations (Huang et al., 2017) and the spatio-temporal modelling of air pollution (Cameletti et al., 2013).
Here our goal was to construct an improved map of terrestrial radiation in Switzerland as a tool for exposure assessment in epidemiological studies. Since the work of Rybach et al. (1997), a more extensive set of radiation measurements has become available. Together with advances in spatial statistics and increased computing power, these should allow the development of more accurate maps of terrestrial radiation than those currently available for Switzerland.

Radiation measurements
We used airborne γ-spectrometry measurements, which were carried out for various purposes including: regular flights to monitor the areas around nuclear facilities (Fig. 2), survey flights to collect reference values in areas of high population density, and training flights for source detection and for international intercomparison exercises. Furthermore, flights traversing large sections of the country as well as targeted flights to observe local anomalies in background radiation have been performed (Fig. 1).
The measurements were performed from a helicopter flying at a height of about 90 m above ground. The system consisted of an a 16.8 litres NaI detector and a spectrometer with 256 channels and energy range of 40-3'000 keV. A spectrum was recorded each second. The ambient dose rate 1 m above ground was computed based on the recorded spectra using the spectrum dose index (SDI) method, described in detail in Bucher (2001). The SDI method has been calibrated by dose rate measurements on ground. The field of view of the detector corresponds to a surface of about 300 × 300 m 2 on ground. This results in a dense coverage and a large overlap of sequential measurements.
The measurements have a relative uncertainty of about 10% for terrestrial radiation. Windows of the spectra allow inferring the soil concentration of 40 K and the elements of the uranium and thorium decay chains, albeit with larger uncertainty. Measurements for 137 Cs mostly lie below the detection limit, however the contribution of 137 Cs to ambient dose rates is incorporated in the estimates for the terrestrial radiation in the SDI method.

Predictors
As potential predictors of ambient dose rates, we considered information on the local geology and land coverage. Geological maps of Switzerland were obtained from the Federal Office of Topography (Swisstopo). We included information on the tectonic plate (19 categories) and on the lithology (5 categories) of surface rock (Federal Office of Topography swisstopo, 2005). Tectonic information captures large scale geological units to which the local geological formations belong. Lithological information refers to the geological processes of rock formation including sedimentary, metamorph or magmatic formation processes as well as categories for unconsolidated rocks and glaciers. Individual categories are listed in tables 4,5 and 6.
Land coverage information was extracted from the areal statistics obtained from the Federal Statistical Office (FSO). Grid cells of 100×100 m 2 are classified into six pre-defined classes based on areal photographs taken during the period of 2004-2008. The six classes refer to artificial surface covers, three vegetation categories (grass, bushes, trees), natural surfaces without vegetation cover and water and wetland surfaces.
We included daily rainfall data (MeteoSwiss, 2013) after the Chernobyl accident as we still expect contributions to the ambient dose rate from 137 Cs contamination from the aftermath of the accident. Most highly affected areas lie in the canton of Ticino. We aggregated the rainfall over the days from 30 April until 5 May 1986. The choice of the days is based on air filter measurements in Fribourg (western Switzerland) and on an animation of the atmospheric spread of the radioactive cloud by the IRSN (Albergel et al., 1988).

Data cleaning and thinning out
We removed values influenced by artificial point sources and measurements at locations where the terrestrial radiation was shielded by water bodies (Fig. 2). Point sources include nuclear power plants, the research site of the Paul-Scherrer Institute (including the intermediate storage facilities for nuclear waste), and a building belonging to the former test reactor facility in Lucens. We investigated spatial correlation with variograms.
Consecutive measurements are correlated due to the overlap of the field of view. This overlap introduces an additional correlation structure on the observation level and masks the inherent correlation of the underlying spatial process (Heersink et al., 2013). We thinned out the measurements by considering only every fifteenth measurement. Details with regard to the thinning out and the reasoning behind it are provided in A.

Model definition
We modelled the ambient dose rate using the following log-linear mixed-effects model where Y (s) is the log-transformed dose rate at location s, β is a vector containing the coefficients of the fixed effects of covariates X(s), U (s) is a Gaussian random field (GRF) with Matérn covariance function and (s) is white noise with variance σ 2 .
The Matérn covariance function has a scale parameter κ > 0 and smoothness parameter ν > 0.
It is defined as where s i − s j is the Euclidean distance between locations s i and s j , K ν the modified Bessel function of the second kind, and σ 2 U is the marginal variance of the spatial process U (s).
Formulated as a hierarchical Bayesian model, we write: where Θ is a vector of parameters of the GRF and Σ its covariance matrix with Σ ij = σ 2 U Cor(U (s i ), U (s j )) and Π denotes the prior distributions.

Inference
Fitting models involving GRFs becomes cumbersome with increasing number of measurements, as the computation involves inverting large matrices (Lasinio et al., 2013). As work-arounds, various approximations have been proposed both within the frequentist and Bayesian framework. These include, for example, covariance tapering (Furrer et al., 2006) fixed-rank kriging (Cressie and Johannesson, 2008), and an approach based on SPDEs (Lindgren et al., 2011). A comparison of proposed techniques suggested similar performance of these methods in a case study competitionHeaton et al. (2019).
We applied the SPDE approach (Lindgren et al., 2011) and fitted the model using integrated nested Laplace approximations (Rue et al., 2009). Integrated nested Laplace approximations (INLA) is a deterministic alternative to Markov chain Monte Carlo (MCMC) methods for fitting (latent Gaussian) Bayesian models (Rue et al., 2009). In the SPDE approach, the Gaussian random field is linked explicitly to a Gaussian Markov random field through the solution U (s) of the SPDE where W is Gaussian white noise, ∆ is the Laplace operator, and U (s) is a continuous GRF with Matérn covariance (Whitle, 1954). Lindgren et al. (2011) represent a weak solution to eq. (4) as a Gaussian Markov random field (GMRF) using the finite element method, expressing the GMRF as linear combination of basis functions defined on the nodes of a triangular mesh. This allows for an explicit link between the GRF and GMRF. A more detailed description can be found in (Krainski et al., 2018). Exact solutions are obtained on the nodes of the mesh and linearly interpolated to a continuous field.
To fit our models, we used the dedicated package R-INLA in the R computing environment (Martins et al., 2013;Lindgren et al., 2015)(www.r-inla.org). In R-INLA, the smoothness parameter ν of the Matérn covariance function is coupled to the parameter choice in the SPDE, and does not need to be additionally set. We chose α = 2, corresponding to ν = 1 (since d = 2). The R-INLA package provides built-in functions to construct the mesh. More nodes, which result in a denser mesh, allow for a smoother representation of the field, but increase the computational burden. The mesh density was tuned by providing the minimal (3.5km) and maximal distances (5km) between nodes together with the minimal angle (31 • ) between edges.
To define priors for the hyper parameters of the spatial field U (s) (σ U and κ = √ 8ν/ρ, where ρ is the (practical) range of the spatial correlation), we used penalized complexity priors (PCpriors) (Simpson et al., 2017). These priors penalize the complexity of the model by shrinking the standard deviation of the field to zero and shrinking the range of the Gaussian field towards infinity. We specified the priors through P (σ U > 10) = 0.01 for the standard deviation of the field, assuming that it is unlikely to have such high spatial variation of radiation left after adjusting for the covariates, and P (ρ > 15) = 0.5 for the range of the field, expecting the range to be in the order of 10-20km after exploring the variograms. Normal priors with mean zero and precision 0.001 were used for the fixed-effects β and an inverse gamma prior with shape 1 and scale 5×10 5 for the precision (1/σ 2 ).

Extended model
To allow for a more complex spatial correlation structure, we extended the model by adding a second spatial random effect U 2 (s). The extended model can be written as where we intended U 1 and U 2 to capture both large scale correlation allowing to make predictions for areas not covered by in areal survey and short range variation helping to fit the data better in densely surveyed areas.

Model comparison
We compared the mixed-effects model and the extended mixed-effects model to simpler models.
We were interested in whether the added complexity of the (extended) mixed-effects model results in a better predictive performance than the simpler model 1. We also considered the following simplified models: • Bayesian linear regression (fitted with INLA): • Spatial Bayesian random-effects model: with Y (s), βX, and U (s) as above. The same set of covariates, parameters for mesh construction, and priors were chosen for the simplified models.
To assess model performance we conducted two types of cross-validation: 1) Random crossvalidation by randomly splitting the data into a training (70 percent) and a validation set (30 percent), and 2) Spatial cross-validation by partitioning the country into grid cells of 15 × 15 km, assigning each cell randomly to one out of four folds and recomputing the model four times, each time leaving out one of the folds. Assignment to the folds was done using the R package blockCV (Valavi et al., 2018).
The spatial cross-validation was tailored to a typical range of extrapolation from surveyed into non-surveyed areas that would be required for exposure assessment in epidemiological studies. The size of the spatial blocks was chosen based on the distance of residential address geocodes to the closest measurement. Geocoded locations of residential buildings in Switzerland were obtained through the Swiss National Cohort (Bopp et al., 2009). We considered locations further than 250 m away from the closest measurement and chose the side length of a block as twice the 90 percent quantile of the distribution of distances between residential locations and closest measurement.
As performance measure for comparing models, we considered the root mean square error defined as where the indices p i and m i stand for predicted and measured values at the locations i = 1, . . . , N for which predictions are made. We considered the mean of the predicted posterior y p as point prediction in order to facilitate the interpretation. The RMSE and the R 2 (defined as Cor(Y p , Y m ) 2 ) are separately reported for the model fit to both the training and validation folds.
In addition, we computed the continuous ranked probability score (CRPS) (Gneiting and Raftery, 2007), a measure of predictive accuracy defined as where

Sensitivity analysis
For the selected model, we assessed whether modifications of the proportion of measurements kept during thinning out (every 10 th and every 20 th vs. every 15 th in the main model) and of the mesh (denser and less dense than main model) affected the results. We compared the resulting predictions on a 1 × 1 km grid by computing the R 2 between the means of pointwise posteriors, and examined changes in the resulting hyper parameters and coefficients.

Predictions
We computed posterior mean, mode, and standard deviation of terrestrial radiation on a 100 × 100 m 2 square grid over Switzerland using the best performing model.
Calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern, and the CX1 Imperial College London cluster.

Description of measurements
The database consisted of 654'530 measurements made during 214 flights. After cleaning, inculding removal of those influenced by point sources or water bodies, 589'197 measurements were available.
The measurements densely cover large parts of the northern part of Switzerland, where population density is highest. The coverage and the partition into 4 spatial folds are illustrated in Figure 1. Variograms before and after including an external trend revealed no (global) directional anisotropy, but indicated spatial processes at different ranges (Fig. 3), motivating the extension of the mixedeffects model with a second spatially correlated GRF.

Comparison of fitted models
The fitted hyper parameters of the different models are shown in Table 1. The practical range was estimated to be 13.9 km in the mixed-effects model, and 15.3 km in the pure spatial model. In the extended model, the fitted practical ranges were 26.6 km for U 1 and 1.7 km for U 2 . The resulting β's are displayed in Tables 4, 5, 6, and 7 (Appendix).
The extended mixed-effects model performed best in both cross-validation settings and across all measures (Table 2). It achieved an R 2 of 0.4 and CRPS of 0.153 in the spatial cross-validation compared to an R 2 of 0.27 and CRPS of 0.217 of the standard mixed-effects model. As expected, all models performed better in a random cross-validation than in a spatial cross-validation setting, where measurements of the validation and training sets are further apart. The differences in performance between the validation settings were more pronounced for the models including spatially correlated random-effects than for the linear model. Moreover, while the linear model is clearly outperformed by the other models in the random cross-validation, only the extended mixed-effects

Sensitivity analysis
When fitting the extended mixed-effects model to different subsets of the data (every 10th, and every 20th measurement included), the resulting predictions did not change substantially (Table 3).
The R 2 between the predictions based on these subsets and the subset used in the main analysis were 0.96 and 0.97 respectively. When fitting the model using a less dense or a denser mesh, we obtained an R 2 larger than 0.99. The fitted hyper parameters are similar (Table 3). Tectonic units that are related to higher dose rates are Tertiaere Intrusiva und Extrusiva (0.43(0.13)), Unterostalpine Decken (0.26(0.08)) and Oberostalpine Decken (0.32(0.10)). Surfaces * CV, Cross-validation; CRPS, continous ranked probability score; RMSE, root mean square error a The data was split randomly into a training (70 percent) and validation (30 percent) sample. b The country was partitioned into 15×15km blocks, which then were assigned randomly to one out of four folds (Fig. 1). The model was recomputed four times, each time leaving out one of the four folds for validation. Displayed are the performance measures averaged over the four folds. covered by bush (-0.019(0.006)) or tree (-0.038(0.004)) vegetation and wetland (-0.13(0.02)) surfaces are associated with lower dose rates, whereas cumulative rainfall after Chernobyl shows a positive association (5.1(1.8)). The contrast between coefficients for the categorical covariates is shrunk in the mixed-effects and extended mixed-effects model compared to the linear model.

Predicted dose rates
Predictions of the preferred extended mixed-effects model, shown in Figure 5, have lower uncertainty in areas not covered by areal survey compared to the simple mixed-effects model (Fig.   4, right panel). Due to the coverage of the biggest cities by measurements, uncertainty in the most densely populated areas tends to be lower than for rural areas. The geographically weighted distribution of predicted values of the preferred model shows a narrow concentration around the value of 50 nSv/h with a heavy tail extending to values over 100 nSv/h (Fig. 6).

Discussion
We provide new estimates of terrestrial gamma radiation for Switzerland based on an spectrometric areal survey of large parts of the territory and information on land coverage, surface rock and underlying geology. The best performing model includes two spatially correlated random effects capturing short and longer range correlations, respectively. This model consistently outperformed a model with only one spatially correlated random effect and the simple linear model, particularly Figure 4: Maps of pointwise posterior mean (predicted mean, left) and standard error of predictions (right) for log-transformed dose rates from terrestrial gamma radiation based on the linear regression, random-effects, mixed-effects and the extended mixed-effects models. Note that the standard error of the predictions refers to the predicted mean surface and does not include the Gaussian white noise (s).
in spatial cross-validation, where training and validation sets were spatially separated.
Areas whose local geology is dominated by sedimentary or unconsolidated rocks generally exhibit lower dose rates than areas with crystalline rocks. The general pattern displayed by all models Figure 5: Backtransformed ambient dose rates from terrestrial radiation predicted by the extended mixed-effects model. Figure 6: Geographically weighted distribution of dose rates from terrestrial gamma radiation predicted by the extended mixed-effects model. is similar to the map of terrestrial radiation developed by Rybach et al. (1997). Contrary to their work, we did not consider the contribution of 137 Cs separately, in view of the levels of radiation stemming from Caesium measured over the last two decades, while the map developed by Rybach et al. (1997) reflects levels measured in 1989 and 1990, i.e. relatively shortly after the Chernobyl accident.
Previous studies in France and Great Britain applied variants of kriging to model the geographic variation of terrestrial radiation. Conceptually our approach is similar, but it is embedded in a Bayesian framework. While we worked with outdoor measurements, both aforementioned studies focused on indoor measurements. Indoor measurements are influenced by building materials and shielding.
Both Warnery et al. (2015) and Chernyavskiy et al. (2016) reported mean square errors (MSE) from a random cross-validation setting of between 360 to 400 (nSv/h) 2 , which is higher than in our work. Kendall et al. (2016) and Chernyavskiy et al. (2016) reported the R 2 in Great Britain to be 0.2-0.27 for different modelling approaches, which is lower than in this work.
Measured dose rates had a similar range but a lower mean in Switzerland (mean: 54.25 nSv/h, range: 7-417 nSv/h) compared to France (mean: 76 nSv/h, range: 13-349 nSv/h) (Warnery et al., 2015), where these were based on dosimeter measurements from veterinary clinics. Consequently, the range of the predicted values is very similar between the two countries. The fitted range of the spatial correlation is roughly one order of magnitude smaller in Switzerland compared to France.
An explanation could be the higher measurement density in our study, allowing us to observe spatial correlations at smaller distances.
A strength of our study is the high resolution of measurements in areas covered by areal surveys.
However, coverage was patchy leaving large non-surveyed gaps. Single measurements observe a field of view of 300 × 300 m 2 and thus measured values need to be interpreted as a weighted averaging of the terrestrial radiation levels over the scale of a few 100 square meters. Variations occurring at a smaller scale cannot be captured.
We did not separately model the contributions stemming from Caesium-137 as most of the measurements were below the detection limit of the measurement device. Decay and vertical migration of 137 Cs since the Chernobyl accident significantly influenced the observed levels of ambient dose rates, most strongly in Ticino. As the measurements were conducted between 1996 and 2018, the coefficient for rainfall is expected to capture an average effect from 137 Cs. As most of the measurements have been performed after 2000 and contamination has been relatively low in the rest of Switzerland, we expected the spatial distribution to be stable enough to neglect the temporal trend without loosing relevant features and predictive performance.
Several aspects should be considered when using the presented maps for exposure assessment in epidemiological studies. Studies for which exposure over the past few decades is of interest should attempt to separately model the temporal trend of the 137 Cs contribution. Furthermore, it must be kept in mind that our map represents ambient doses outdoors. While this can serve as a proxy measurement for indoor exposure, further research on the differences between indoor and outdoor exposure levels in Switzerland might allow improved estimation of indoor exposure. Due to the gaps in the data coverage, there are areas of relatively large uncertainty with regard to the predicted dose rates. In addition to the map, the chosen modelling approach also allows producing maps of the predictive uncertainty. These could be used to propagate the uncertainty forward to models linking exposure with disease outcomes so as to be correctly reflected in estimates of potential health effects. Besides the exposure assessment in epidemiological studies, the provided map can serve as a baseline map for assessing potential contamination.
Our finding that two spatial components with differing ranges improved prediction suggests that spatial variation of terrestrial radiation occurs at different scales. The spatial correlation of terrestrial γ-radiation might thus be better captured by multiple processes acting over different ranges than by a single spatial process. Models that are able to incorporate such behaviour might prove most suitable to map the variations of terrestrial radiation. This also offers the possibility of improving models in future. Similar situations may exist for other exposures.
The map, formatted as shapefile, and the R code used for estimating the model are available on Github (https://github.com/FollyCh/TerrMapCH). Also we provide 100 realisations drawn from the joint posterior distribution that can be used for error propagation. and assisting our selection of covariates.

Declaration of competing interest A Thinning out
Our decision to thin out the measurements used for model fitting were based on the following considerations.
1. Overlapping observation windows: consecutive measurements cover partly overlapping areas.
Additionally, the measurement errors might be correlated between subsequent measurements, which would introduce spurious correlations.
2. Extreme extrapolations: when we fitted the models to the full data set, we observed extreme values for extrapolations into unobserved areas, most pronounced around measured local peaks. We assume the reason for these extremes to be trends at the borders of the surveyed areas propagated some distance into the areas not surveyed.

Redundant information:
As the measurements densely cover surveyed areas, subsets of the data might still contain all relevant information. By thinning out we can save computation time and disk space. Aggregating or thinning out should lead to more stable computations and potential impacts of spurious correlations are mitigated.
Our decision to select only every 15th measurement was based on variograms of consecutive measurements in the directions of flight paths. To compute these, we defined the x-coordinate as the enumeration of subsequent measurement points and the y-coordinate as zero and calculated separate variograms for randomly chosen flights. Results indicated a strong serial correlation of measurements along flight paths which approaches zero only after a lag of 15 to 20 subsequent measurements.      Table 7: Fitted Coefficients in the different modelling approaches for cumulative rainfall after Chernobyl. We rescaled the cumulative rainfall to meters to compute the models.