*Article* **From Spatial Characterisation to Prediction Maps of the Naturally Occurring Radioactivity in Groundwaters Intended for Human Consumption of Duero Basin, Castilla y León (Spain)**

**David Borrego-Alonso 1,\* , Antonio M. Martínez-Graña <sup>2</sup> , Begoña Quintana <sup>1</sup> and Juan Carlos Lozano <sup>1</sup>**


**Abstract:** According to the European Council Directive 51/2013 EURATOM, the radionuclide content in human consumption waters must comply with the stated recommendations to ensure the protection of public health. The radiological characterisation assessed in Laboratorio de Radiaciones Ionizantes y Datación of Universidad de Salamanca, in more than 400 groundwater samples gathered from intakes intended for human consumption from the Castilla y León region (Spain), has provided worthwhile data for evaluating the spatial distribution of the radioactivity content in the Duero basin. For this purpose, geostatistical exploration and interpolation analysis, using the inverse distance weighting (IDW) method, was performed. The IDW prediction maps showed higher radioactivity occurrence in western and southern areas of the study region, mainly related to the mineralogical influence of the igneous and metamorphosed outcroppings in the Cenozoic sediments that formed the porous detritic aquifers of the Duero basin edges. The hydraulic characteristics promote the distribution of these radioactivity levels towards the centre of the basin as a consequence of the unconfined nature of the aquifers. Prediction maps provide a worthwhile tool that can be used for better planning and managing of groundwater monitoring programmes.

**Keywords:** natural radioactivity; spatial distribution; IDW; prediction maps; groundwater; drinking water

## **1. Introduction**

It is widely known that the quality assessment of groundwater arouses great interest in environmental management for planning the use of water. Groundwater is not only a strategic resource for drinking water supply for human consumption, but also to satisfies irrigation demand, especially in arid and semi-arid regions, where intensive agriculture is one of the main economical driving forces and the seasonal variability in the precipitation regimes frequently lead to drought periods. In Spain, groundwater withdrawal is estimated at 22% of the total [1], However, there exists a wide variability range among the different administrative authorities commissioned for hydrographical planning. In Castilla y León the primary groundwater abstraction comes from the Duero Hydrographical Confederation (DHC), which extends through 82.2% of their territory, and is approximately 1220 hm3/yr, mainly intended for agriculture and livestock farming (78%) and domestic human consumption (16%). Their use becomes especially significant in regions where the superficial supply can no longer provide enough water, due to the absence of reservoirs or during the aforementioned seasonal droughts. Hence, it is essential to accomplish exhaustive characterization assessments that ensure the quality levels of the groundwater resources intended for human consumption are in compliance with public health recommendations.

**Citation:** Borrego-Alonso, D.; Martínez-Graña, A.M.; Quintana, B.; Lozano, J.C. From Spatial Characterisation to Prediction Maps of the Naturally Occurring Radioactivity in Groundwaters Intended for Human Consumption of Duero Basin, Castilla y León (Spain). *Agronomy* **2022**, *12*, 2059. https:// doi.org/10.3390/agronomy12092059

Academic Editor: Paul A. Ty Ferré

Received: 26 July 2022 Accepted: 26 August 2022 Published: 29 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Given their critical impact on human health, the determination and spatial distribution of heavy metal concentrations have been addressed in different environments [2–5]. The naturally occurring concentration of these elements in groundwater depends on the involved processes in the transference from rock and soil minerals to aquifers. Apart from the chemical toxicity associated with the ingestion of heavy metals through the food chain and drinking water, some of them, such as uranium, radium, lead or polonium, are also radiotoxic due to the ionising radiation emitted as a consequence of their radioisotope disintegration. Consequently, regulatory frameworks developed during recent years, have increasingly focused on the radioactivity content in groundwaters intended for human consumption to protect public health from radionuclide substances.

The naturally occurring radionuclides present in groundwaters come from the three natural decay chains that form part of the crystalline structure of the bedrock minerals. Taking into account their relevance as a key resource for water availability in many worldwide regions, several works [6–10] have focused on the hydrogeochemical processes involved in the radionuclide mobilisation from mineral-bearing to the groundwater environment. Background radioactivity levels on groundwaters are often innocuous, and it is also a fact that there are regions prone to lithological features with specific hydrogeochemical settings where the radionuclide presence is enhanced [11–15], highlighting the need for quantifying and monitoring of their content. Prolonged radiation exposure due to the ingestion of drinking water containing certain radionuclide levels is likely to lead to an increase in the risk of stochastic effects on the development of severe diseases [16], related to the radiological toxicity and the chemical behaviour of uranium, radium, polonium and lead. Hence, an accurate radiological characterisation of the groundwater bodies is an important challenge to guaranteeing safe human water consumption. The current legal framework for the protection of public health regarding the radioactivity parameters was laid down in the European Council Directive 51/2013/EURATOM [17]. Spain, as a Member State of the EU, also incorporates the European Directive requirements in its legislation to ensure compliance with the quality radioactivity standards [18]. According to the safety guidelines of radioactivity in drinking water, the control of the indicative dose (*ID*) must be addressed by determining the activity concentrations of 226,228Ra, 234,238U, <sup>210</sup>Pb and <sup>210</sup>Po, as established in Annexe X of the Royal Decree (RD) 314/2016. Once the water has been radiologically characterised, the gross α- and gross β-activity determinations can be used as a screening strategy for monitoring the radioactivity levels. If gross α- and gross β-activity values are kept below their screening levels, 0.1 and 1.0 Bq/L, respectively, it may be assumed that the *ID* does not exceed the parametric value of 0.1 mSv, and no further costly and expertise-specific radiological determinations are necessary. Between 2017 and 2021, in the framework of the project promoted by the General Directorate of Public Health of Castilla y León, our laboratory carried out a comprehensive radioactivity characterisation of the groundwater intakes placed throughout the whole region, which comprised an analysis of the 238,234,235U, 228,226Ra, <sup>210</sup>Pb and <sup>210</sup>Po activity concentrations for the *ID* assessment. Furthermore, their respective gross α- and gross β-activities, hereinafter referred to as *a*<sup>α</sup> and *a*β, respectively, were also determined, due to their wide use in water radioactivity monitoring.

In the last decade, geostatistic and geospatial analysis has become a highly useful and versatile tool for understanding environmental variability and the spatial distribution of groundwater quality [19–22]. Some publications have reported that high radioactivity levels on groundwaters are severely associated with the bedrock type and the hydrogeochemical variables [22–25]. Geographic Information System (GIS) and geostatistical tools implemented in the software enable the assessment of spatial potential relationships between the groundwater radionuclide content and the environmental setting [26,27]. Despite the existence of data-driven methods used for environmental studies, such as artificial neural networks [28], we performed the inverse distance weighted (IDW) interpolation method, which allows the estimation of the unmonitored area based on the quantitative measurements carried out in neighbourhood sampled locations, to assess the spatial distribution of the radioactivity content in the aquifers of Castilla y León. The IDW interpolation method is widely used in geospatial studies related to groundwater pollution and quality assessment [29]. Despite the different interpolation methods, the most appropriate depends on the type of data parameters and their predictability is determined by multiple factors related to the distribution of the sampling points [30]. Some publications have addressed the monitoring, geospatial analysis and projecting of prediction maps of the environmental radioactivity in waters, such as <sup>222</sup>Rn in southern Belgium [31] or *a*<sup>α</sup> and *a*<sup>β</sup> in Turkey [32]. Our study presents a radioactivity assessment in a region where different main lithologies and permeabilities define the aquifers of Castilla y León. Given that the prediction maps of the radioactivity parameters in the study region constitute an important milestone, not only for a better understanding of the groundwater radionuclide distribution, but also because they provide a powerful tool from a management perspective for compliance with the current law. Thus, the current study was conducted to assess the spatial distribution of the main naturally occurring radionuclides using the IDW method, by considering analysed samples from hydraulically interconnected aquifers belonging to the DHC, and to provide groundwater prediction maps that can be used as a tool in the groundwater planning.

## **2. Materials and Methods**

#### *2.1. Study Area*

Castilla y León (CyL) is situated in the NW of the Iberian Peninsula, and occupies approximately 94,225 km<sup>2</sup> , thus, constituting the largest region of Spain and the third largest of the European Union.

Geologically, regarding the lithostratigraphic and structural features, different domains can be distinguished, as can be observed in Figure 1. The Hercynian basement is comprised mainly of igneous and metamorphic outcroppings, originating during Precambrian and Paleozoic, and constitutes the Central System (CS). Throughout the westernmost areas of the provinces of Zamora and Salamanca and the southern areas of Salamanca, Ávila, Segovia and León, acidic plutonic outcroppings and Precambrian-Cambrian metasediments, which mainly consist of slates, greywackes and sandstones, and small occurrences of quartzites, conglomerates and black-carbonaceous slates, mainly comprise the basement, which constitutes the Schist-Greywacke Complex (SGC). Igneous intrusive bodies are broadly formed by two-mica granites and other biotitic granitoids. Throughout the western Asturian-Leonese zone, and in small outcrops located in the Iberian System, located in the eastern region of CyL, Cambrian-Ordovician-Silurian sediments, with some interlaced levels of carbonates, overlying in discordance with the Precambrian metasediments, such as slates, shales, quartzites and *Ollo de Sapo* gneiss rocks predominate. The Cantabrian zone, which occupies the northernmost areas of the provinces of León and Palencia, is comprised of Precambrian and Paleozoic sediments deposited in shallow marine environments, originating in limestones and sandstones. Terrigenous and carbonated Mesozoic materials mainly constitute the basement of the Basque-Cantabrian Basin, which extends throughout the north of Palencia and Burgos provinces, and the Castillian part of the Iberian Range, represented mainly by the Cameros Basin between SE of Burgos and the north of Soria. The sedimentary sequence of the different sectors of the Cenozoic Duero Basin (DB) is influenced by the bedrock composition that surrounds the Spanish Northern Plateau. It is worth pointing out the significant impact on the composition of the siliciclastic rocks, that fill the southern and western sectors of the DB, of acidic intrusive rock debris and metamorphic sediments belonging to the CS. Finally, fine to coarse detrital grains (sands, clays, limes and conglomerates), associated with the current riverine systems, constitute the Quaternary deposits, which originated from the drainage of the DB materials and the continuous denudation of Mesozoic mountain edges [33].

denudation of Mesozoic mountain edges [33].

belonging to the CS. Finally, fine to coarse detrital grains (sands, clays, limes and conglomerates), associated with the current riverine systems, constitute the Quaternary deposits, which originated from the drainage of the DB materials and the continuous

**Figure 1.** Map of the main lithologies and sampling sites (in black dots) of the studied groundwaters of Castilla y León (Spain), belonging to the centre of DHC. **Figure 1.** Map of the main lithologies and sampling sites (in black dots) of the studied groundwaters of Castilla y León (Spain), belonging to the centre of DHC.

The spatial distribution of the natural radionuclides, the groundwater occurrence of which is controlled by the processes within the water–rock interactions, is, thus, highly influenced by the main lithologies represented by granitic and metamorphic outcroppings, detritic materials and carbonate rocks and their respective abilities to store, transmit and

discharge groundwater [34,35]. Accordingly, 20 hydrogeological unities are defined in the centre of the DHC which are formed by groundwater masses or aquifer systems that present similar hydraulic, hydrogeological and distribution characteristics and also share the same environmental planning goals, as illustrated in Figure 2. transmit and discharge groundwater [34,35]. Accordingly, 20 hydrogeological unities are defined in the centre of the DHC which are formed by groundwater masses or aquifer systems that present similar hydraulic, hydrogeological and distribution characteristics and also share the same environmental planning goals, as illustrated in Figure 2.

The spatial distribution of the natural radionuclides, the groundwater occurrence of which is controlled by the processes within the water–rock interactions, is, thus, highly influenced by the main lithologies represented by granitic and metamorphic outcroppings, detritic materials and carbonate rocks and their respective abilities to store,

*Agronomy* **2022**, *12*, 2059 5 of 19

**Figure 2.** Permeability map of the Castilla y León region according to the main lithologies. Black lines limit the 20 hydrogeological units of the DHC studied. **Figure 2.** Permeability map of the Castilla y León region according to the main lithologies. Black lines limit the 20 hydrogeological units of the DHC studied.

Regarding permeability, there are broadly two types of productive aquifer systems associated with the rock type and permeability characteristics. On the one hand, the intergranular porous aquifers are formed in detritic materials, such as sands, gravels and conglomerates, and extend throughout the Cenozoic sediments of the DB and tectonic depressions of Ciudad Rodrigo or Amblés. These aquifers are included in the following Regarding permeability, there are broadly two types of productive aquifer systems associated with the rock type and permeability characteristics. On the one hand, the intergranular porous aquifers are formed in detritic materials, such as sands, gravels and conglomerates, and extend throughout the Cenozoic sediments of the DB and tectonic depressions of Ciudad Rodrigo or Amblés. These aquifers are included in the following hydrogeological unities: Obrigo, Esla-Cea and Cea-Carrión "Rañas", Esla-Valderaduey Region, Páramo de Torozos, Duero Central Region, Burgos-Aranda, Alluvials of Duero River,

Almazán Depression, Arenales, Ciudad Rodrigo-Salamanca, Valdecorneja and Amblés Valley. Theyepresent 85% of the territory and 70% of the total groundwater volume of the most important hydric resources. The Cenozoic aquifer systems are formed by a wide sands lens with a semipermeable lime-clayish matrix surrounding them. These aquifers behave as a great anisotropic interconnected aquifer complex, which presents a variable thickness of up to 2000 m depth in the basin centre. Infiltration and sideway flow are the main recharge mechanisms. Generally, the groundwater flows towards the basin centre through the current riverine drainage system. Productivity and transmitivity are highly variable, 10–40 L/s and 5–100 m2/d, respectively. The superficial (5–15 m), low-transmissive and low-productive porous aquifers formed in the Quaternary sediments, such as rañas, alluvials and sandbanks, are, likewise, variable. On the other hand, the fissured carbonated aquifers, formed in the Mesozoic and Cenozoic limestones and dolomites, present a permeability related to fissures due to dissolution and structural processes. Mesozoic formations, which comprise the northern and western peripheral areas of the DB, constitute broadly unconfined karst aquifers recharged by freshwater runoff and interconnected groundwater lateral flow with other units, and their transmitivity and productivity range between 20–1500 m2/d and 5–140 L/s, respectively. In the centre of the DB, limestone formations extend, giving rise to the Páramos (Cúellar, Duratón and Torozos) that constitute unconfined aquifers with a subhorizontal tabular structure, 5–50 m depth, radial flow and recharge– discharge through several springs, as well as infiltration. It can be considered that to a greater or lesser extent, all these aquifers are interconnected. It must be noted that the igneous and metasedimentary outcroppings, which extend throughout the Central System and the North-Cantabrian zone, give rise to limited and fractured regionally insignificant aquifers, hydraulically confined, and characterised by very low permeability. In the Ebro Basin, which occupies the Basque-Cantabrian zone, Cenozoic basin sediments exist with clayey interlaced levels that limit the hydraulic relation but give rise to some confined aquifers or locally significant aquifers.

### *2.2. Groundwater Samples*

Between January 2017 and November 2021, two hundred and forty-four drinking water samples were gathered from groundwater intakes situated within the hydraulically interconnected DHC. The georeferenced location of the sampling points is shown in Figure 1. The selection of the sampling points responded to strategic goals stated by the Dirección General de Salud Pública (DGSP) of Junta de Castilla y León to comply with the guidelines for the radiological characterisation of the groundwater masses intended for human consumption.

The water samples were mainly collected by the Pharmaceutical official service of the DGSP in three 10-L PTE vessels and then transported to Laboratorio de Radiaciones Ionizantes y Datación of Universidad de Salamanca for measurement and analysis to determine the radionuclide content. Immediately after this measurement, water samples were acidified with 1 mL/L of analytical grade 65% HNO<sup>3</sup> for proper preservation.

## *2.3. Radioactivity Measurements*

The activity concentration of 238,235,234U, 228,226,224Ra, <sup>210</sup>Pb and <sup>210</sup>Po were determined using γ-ray and α-particle spectrometry techniques. The *a*<sup>α</sup> and *a*<sup>β</sup> values were determined using the optimised thin source deposit and proportional counting (TSD-PC) method [36].

The γ-sources were prepared by the thermal concentration of a twenty-five-litre aliquot up to 50 mL final volume. The concentrated samples were then transferred to cylindrical 52 mm diameter and 48 mm height PTE beakers. Two low-level background HPGe detectors, manufactured by Canberra, BE5030 and GR2520 models, hereinafter referred to as BEGe and REGe, respectively, were used for the γ-ray spectrometry measurements. Both detectors were continuously ventilated with N2. The BEGe detector was made of a p-type high-purity Ge crystal of broad-energy range with an active volume of 117 cm<sup>3</sup> . Its relative energy peak efficiency at 1332 keV was 50% and its nominal resolution at 122 and 1332 keV

were 0.75 and 2.20 keV, respectively. The detector was surrounded by a passive shielding consisting of 10 cm thick iron melted and 5 cm thick old lead, which was internally lined with 2 mm thick high-purity electrolytic Cu. The REGe detector was made of an n-type high-purity Ge crystal with an active volume of 59 cm<sup>3</sup> . Its relative peak efficiency at 1332 keV was 22.60% and nominal resolution at 122 and 1332 keV were 0.842 and 1.79 keV, respectively. The detector was surrounded by a passive shielding made of 27.5 cm thick iron melted, internally lined with 2 mm thick high-purity electrolytic Cu. The data acquisition was done with coupled electronics consisting of an integrated module Canberra DSA1000 model, including a 16 K multichannel analyser. The γ-sources were measured for times ranging between 250,000 to 450,000 seconds. For efficiency calibration, both Monte Carlo (MC) simulation and experimental methods were used [37,38].

The acquired γ-spectra were analysed using the Gamma-Live Expert Analyser, GALEA, developed in our laboratory mainly for the analysis of spectra shaped by the natural radionuclide emissions [39]. It included the algorithm COSPAJ for the full continuum spectra fitting [40]. A genetic algorithm was also included that provided the best peak fit, even in multiplets with more than three peaks. The γ-ray emission energies and probabilities were taken from the Nuclide-Lara Library [41].

The 238,235,234U activity concentrations were also determined using α-particle spectrometry. Sources were prepared through the radiochemical separation procedure which involves sequential extraction using UTEVA resin (Triskem) of purified uranium eluate and electrodeposition onto stainless steel discs at 1.8 A for 1 hour [42,43]. The method used for <sup>210</sup>Po determination was based on co-precipitated with Fe(OH)<sup>3</sup> and auto-deposition onto silver discs for 4 hours [44,45]. The activity concentration of <sup>210</sup>Po corresponded to the <sup>210</sup>Po in excess, thus, the activity concentration measured from α-particle spectrometry was corrected with the γ-ray spectrometry <sup>210</sup>Pb activity. Measurements were performed using a PIPS semiconductor detector of 450 mm<sup>2</sup> active area, Ortec BR-SNA-450-100 models, housed in an Ensemble Spectrometer with Alpha-Duo Modules (Ortec) for uranium sources and Canberra A450-18 AM model, coupled to low-noise preamplifiers and amplifiers, all of them housed in width NIM spectrometers Canberra 7401VR model, coupled to a multichannel analyzer Ortec 920E EtherNIM 16-Input Multichannel Buffer model, was used for polonium. Spectra and hardware were managed by the MCA Maestro Emulator and the suite Alpha-Vision both by Ortec.

According to the European framework, the water quality intended for human consumption should ensure safe radioactivity exposure for public health. Thus, the indicative dose (*ID*) was determined in all groundwater samples, following Equation (1):

$$ID = \sum\_{i=1}^{n} \frac{c\_i(cal)}{c\_i(der)} \le 0.1 \text{ mSv} \tag{1}$$

*ci (cal)* is the activity concentration determined by a radionuclide (*i*) present in the sample, *ci (der)* is the derived concentration of a radionuclide (*i*) listed in RD 314/2016, and *n* is the number of radionuclides with concentration activities higher than the minimum detectable activity (MDA). The contribution of <sup>210</sup>Pb, <sup>210</sup>Po, <sup>226</sup>Ra, <sup>228</sup>Ra, <sup>238</sup>U and <sup>234</sup>U were considered for the *ID* assessment. The value of the derived concentration is based on the radiotoxicological properties of the radionuclides [46].

The *a*<sup>α</sup> and *a*<sup>β</sup> values were ascertained by following the optimised TSD-PC method developed in our laboratory for drinking water [36]. For each sample, efficiencies were directly calculated by spiking with natural uranium standard solution. Measurements were performed in a low-background gas-proportional counter (model LB770, Berthold Technologies). This equipment enables the simultaneous measurement of the α and β emissions. Underground location, the anticoincidence guard detector and the passive shielding based on thick lead blocks ensured low background.

#### *2.4. Statistical Analysis of the Samples*

The spatial and geostatistical analysis was performed using ArcGIS software (v 10.9) which included the 3D Analyst Toolkit for applying the interpolation methods and the semivariogram analysis. The inverse distance weighted (IDW) interpolation technique determined the unsampled locations from a linear weighted function based on the inverse distance between the sampled location, in which values were known. This interpolation method takes into account that the nearest measured values have more weight in the prediction than those further away. The predicted surface area was computed using the Equation (2):

$$Z\left(\mathbf{S}\_{\bullet}\right) = \sum\_{i=1}^{N} \lambda\_{i} Z\left(\mathbf{S}\_{i}\right),\tag{2}$$

where *Z* (*Si*) is the measured value in the sample point *i*; *λ<sup>i</sup>* is the unknown weighting for the measured value in the sample point *i*; (*SO*) is the location of the prediction; and *N* is the number of measured values [47].

The weight of the set of measured points in the prediction model diminishes with the distance. Despite other non-deterministic interpolation techniques, such as geostatistical techniques like Kriging, often being applied, the current work considered the IDW technique to be most proper interpolation technique to also consider the spatial distribution of the sampled locations and build prediction surface maps with some measure of certainty and predictive accuracy [48]. The absence of a higher density of the grid sampled points in such an extensive territory led to discarding of Kriging interpolation.

The IDW interpolation method was conducted within the boundaries of the hydraulically connected aquifers. It was adjusted to 12 site points with a variable search radius, and with a grid size of 50 m. The power parameter, which determines the significance of the model, was defined with a value of 2. The break values were established to represent in red colour the interpolated areas where the activity concentration of the specific radionuclide exceeded their contribution to the *ID* parametric value stated. Despite <sup>210</sup>Po also being considered for *ID* calculation, we dide not apply IDW interpolation to it because there were a lot of water samples having <sup>210</sup>Po activity concentration below MDA. Moreover, only 9 samples exceeded 10 mBq/L, which is a magnitude order lower than the activity concentration that leads to presenting an *ID* higher than the stated parametric value, Table 1. In the case of the *a*<sup>α</sup> and *a*β, classification shown in the map legends was performed regarding whether values exceeded the screening value or not, established at 0.1 Bq/L and 1.0 Bq/L, respectively [17,18].

**Table 1.** Derived concentrations corresponding to the main natural radionuclides present in the studied samples for *ID* calculation according to Annexe X of RD 314/2016.


#### **3. Results**

Exploration data analysis was performed to gain a better understanding of the spatial distribution, as well as to identify outliers and trends. For each radionuclide and parameter studied, the histogram provided all the statistics and the frequency distribution shape. The right tail shown in Figure 3a indicated a relative lack of large activity concentration values, evidencing that they were not normally distributed. All of them followed a log-normal decreasing continuous probability distribution as can be seen in Figure 3b, where the QQ plot for *a*<sup>α</sup> is shown, in which normality was then demonstrated, given that most of the data were close to the 45-degree line. In Table 2, where non-log-normal transformed statistic values are summarised, the skewness values, as well as the differences observed between the mean and median, indicated high asymmetry of the data distribution and the interpolation method improved the modelisation of the dataset A three-dimensional view of the data was given by the trend analysis tool which allowed us to identify global spatial

trends in the dataset. According to the distribution pattern observed in Figure 3c, the spatial autocorrelation of high values increased longitudinally from east to west and latitudinally southward, but the latitudinal component had greater weight than the longitudinal one. The variability of the spatial autocorrelation of the sampling points was analysed through the semivariogram cloud (Figure 3d). The theoretical base was the first Tobler's law, which states that the nearest points are more similar than those that are outermost and are then more predictable and less variable. On the *x*-axis the distance separating each pair of points was plotted, while the *y*-axis refers to the relative semivariogram values, which corresponded to the squared difference between each pair of sampling points. As the distance between the pair of sampling points increased, the semivariogram value also increased. The semivariogram cloud shown in Figure 3d indicated that most of the data were spatially autocorrelated. **Table 2.** Statistical analysis of activity concentration of the main natural radionuclides and radioactivity parameters in the interconnected groundwaters of Castilla y León. **234U 238U 226Ra 228Ra 210Pb** *ID a***<sup>α</sup>** *a*<sup>β</sup> Mean (̅) 131 58.2 24.04 9.62 16.1 20.8 215 228 Median (Me) 51.4 21.8 5.30 4.10 4.84 11.6 108.0 132 Minimum (min) 0.07 0.04 0.13 0.03 0.37 0.33 0.67 1.83 Maximum (max) 1616 705 1162 196 512 318 2407 1547 Standard deviation (S) 208 91.1 86.80 18.9 47.4 28.5 323.3 254 Kurtosis 20.0 16.8 101.08 58.4 64.7 53.1 20.184 9.45 Skewness 3.57 3.25 9.16 6.59 7.14 5.72 3.64 2.35

Figure 3d indicated that most of the data were spatially autocorrelated.

given that most of the data were close to the 45-degree line. In Table 2, where non-lognormal transformed statistic values are summarised, the skewness values, as well as the differences observed between the mean and median, indicated high asymmetry of the data distribution and the interpolation method improved the modelisation of the dataset A three-dimensional view of the data was given by the trend analysis tool which allowed us to identify global spatial trends in the dataset. According to the distribution pattern observed in Figure 3c, the spatial autocorrelation of high values increased longitudinally from east to west and latitudinally southward, but the latitudinal component had greater weight than the longitudinal one. The variability of the spatial autocorrelation of the sampling points was analysed through the semivariogram cloud (Figure 3d). The theoretical base was the first Tobler's law, which states that the nearest points are more similar than those that are outermost and are then more predictable and less variable. On the *x*-axis the distance separating each pair of points was plotted, while the *y*-axis refers to the relative semivariogram values, which corresponded to the squared difference between each pair of sampling points. As the distance between the pair of sampling points increased, the semivariogram value also increased. The semivariogram cloud shown in

*Agronomy* **2022**, *12*, 2059 9 of 19

**Figure 3.** Geostatistical exploration corresponding to the *a***α** analysis: (**a**) histogram; (**b**) logtransformed QQ normal plot; (**c**) trend analysis; (**d**) semivariogram. **Figure 3.** Geostatistical exploration corresponding to the *a*<sup>α</sup> analysis: (**a**) histogram; (**b**) log-transformed QQ normal plot; (**c**) trend analysis; (**d**) semivariogram.

**Table 2.** Statistical analysis of activity concentration of the main natural radionuclides and radioactivity parameters in the interconnected groundwaters of Castilla y León.


### *3.1. Radionuclide Analysis and Their Relationship with the Lithological Context and Permeability*

IDW provided the prediction maps for the main natural radionuclides as well as the *a*α, *a*<sup>β</sup> and *ID*, which are represented in Figure 4a–f. The spatial distribution pattern predicted for uranium radionuclides, <sup>234</sup>U and <sup>238</sup>U, showed a similar defined area, as can be seen in Figure 4a,b, but with a greater occurrence of <sup>234</sup>U. Although in the mineral grains both radionuclides were frequently in secular equilibrium, the presence of <sup>234</sup>U in groundwater was often much higher than the <sup>238</sup>U concentration because of the different physicochemical processes involved in their mobilisation. Given that uranium solubility is enhanced under oxidising conditions at near-neutral pH, typical in shallow aquifers, the weathering of mineral-bearing host rock releases the uranium nuclides to groundwater, which is the only mechanism controlling <sup>238</sup>U release from crystal lattice to water phase. However, when <sup>238</sup>U decays by alpha emission the recoil process may eject <sup>234</sup>Th into the mineral lattices close to the water–rock interaction boundary, producing recoil damage tracks where is located. The rapid decay of <sup>234</sup>U from <sup>234</sup>Th (half-life 24.1 d), promotes the <sup>234</sup>U preferential leaching [49].

As seen in Figure 4a, <sup>238</sup>U activity concentrations were higher than 100 mBq/L, which corresponded to the recommended screening level for *a*α, and had been predicted in the southern and central aquifers of the DB region. The concentrations of <sup>238</sup>U extended mainly throughout the hydrogeological units of Los Arenales and Segovia, from the southern edges towards the centre of the DB, reaching values up to and exceeding 700 mBq/L in some points (Table 3). Furthermore, in the surrounding areas of the Duero river, the presence of <sup>238</sup>U was predicted in significantly high concentrations, mainly in the south of the carbonate aquifers of Páramo de Torozos and Duero Central region, as well as in the detritic hydrogeological units of Alluvial of the Duero and Esla-Valderaduey.




**Table 3.** *Cont.*

Similarly to <sup>238</sup>U, the spatial distribution of high <sup>234</sup>U concentrations was found in the southcentral aquifers of the DB, mainly belonging to Los Arenales, Esla Valderaduey, Ciudad Rodrigo-Salamanca and throughout the course of the Duero river, as illustrated in Figure 4b. <sup>234</sup>U concentrations exceeding 100 mBq/L were mainly associated with detritic aquifers, but also some groundwaters hosted in the carbonate rocks developed in the Páramo de Torozos, Páramo de Duratón and Páramo de Cúellar, and the outcroppings in the central-north area of Duero Central region.

The outcropping of U-rich granitoid and other intrusive igneous rocks, as well as metamorphosed Hercynian rocks throughout the CS, gave rise to limited and confined fractured aquifers with <sup>238</sup>U and <sup>234</sup>U concentrations that could even exceed 1300 and 1900 mBq/L, respectively, in some points of these confined aquifers, from previous analysis performed in our laboratory. High naturally-occurring uranium concentrations in soils developed on the Hercynian basement have been reported in the west of Salamanca province [2]. The freshwater infiltration through the fractured basement, runoff in the mountain system of Sª de Ávila and Gredos, and streamflow infiltration from surface bodies recharge the local confined aquifers, where oxic conditions promote the uranium leaching from mineral-bearing to groundwater, leading to an increase in the uranium concentration. The uranium spatial distribution predicted throughout these aquifers could indicate that the Cenozoic detritic deposits comprising the southern part of the DB are influenced by the U-rich mineral composition characteristic of the CS. Porous and moderate to highly productive, predominantly oxic aquifers [50], Cenozoic and Quaternary aquifers of the DB, and the Mesozoic fissured and moderately productive carbonate ones, are recharged through surface infiltration of runoff freshwater, which easily promotes the leaching and dissolution of U-minerals. Furthermore, the general groundwater flux from recharge mountain areas towards the Duero River, as well as the streamflow of fluvial drainage across the Quaternary deposits, may also positively affect uranium mobilisation to groundwater. According to these hydrogeological settings, the IDW method provided a reliable prediction map of the uranium concentrations in the DB aquifers.

**Figure 4.** Spatial distribution of the main natural radionuclides in groundwater and the *ID*. IDW interpolation technique was applied: (**a**) 238U activity concentration; (**b**) 234U activity concentration; (**c**) 226Ra activity concentration; (**d**) 228Ra activity concentration; (**e**) 210Pb activity concentration; (**f**) **Figure 4.** Spatial distribution of the main natural radionuclides in groundwater and the *ID*. IDW interpolation technique was applied: (**a**) <sup>238</sup>U activity concentration; (**b**) <sup>234</sup>U activity concentration; (**c**) <sup>226</sup>Ra activity concentration; (**d**) <sup>228</sup>Ra activity concentration; (**e**) <sup>210</sup>Pb activity concentration; (**f**) *ID*. The hydrogeological unit names corresponding to the numbers illustrated are referred to in Figure 2.

Despite <sup>226</sup>Ra (half-life 1.6 <sup>×</sup> <sup>10</sup><sup>3</sup> y) and <sup>228</sup>Ra (half-life 5.75 y) having the same chemical behaviour, they come from different disintegration chains, <sup>238</sup>U and <sup>232</sup>Th decay series, respectively. Thus, the head chain radionuclide concentrations in aquifer rocks and their widely varying half-lives can affect the different occurrences of <sup>226</sup>Ra and <sup>228</sup>Ra in groundwater [51,52]. As seen in Figure 4c, the interpolation map of <sup>226</sup>Ra occurrence predicted groundwater concentrations between 100 to 500 mBq/L, values which correspond to *a*<sup>α</sup> screening level and the maximum reference concentration (MRC) in drinking waters, respectively, extending through three main locations:


According to the mechanism that may control the presence of <sup>226</sup>Ra in groundwater, the IDW method provided a consistent spatial distribution. The diverse lithological context and the occurrence of high <sup>226</sup>Ra concentration may be mainly due to both mineralogical bedrock composition and the chemical properties that govern the water–rock interaction. On the one hand, the elevated uranium content observed in felsic rock and shale groundwater masses enhances the entry of <sup>226</sup>Ra into the water by the decay of dissolved parent radionuclides, alpha recoil effect, and desorption from mineral grains. On the other hand, the analogue chemical behaviour of radium with other crustal cations, such as Ca or Ba, may lead to it being included in the carbonate complexes and then mobilised to water through dissolution, adsorption or ion exchange with the reactive surface of other colloids present in groundwater. It can be observed in Figure 4c that in the moderately productive carbonate aquifers, in particular, which extended through the Duero Central region, Páramo de Cúellar and Páramo de Duratón, the <sup>226</sup>Ra concentration predicted was <25 mBq/L.

Regarding the <sup>228</sup>Ra concentration, given that it decays from <sup>232</sup>Th, which in most aquifer conditions tends to precipitate as insoluble hydroxides, high values exceeding MRC (200 mBq/L), were not expected in the study area. As shown in Figure 4d, the prediction map showed a co-occurrence of significant <sup>228</sup>Ra concentration in the same areas where <sup>226</sup>Ra-rich was found, corresponding to the southwest of the Esla-Valderaduey region aquifers and in the western of the Ciudad Rodrigo-Salamanca unit, with maximum concentrations of 196 and 164 mBq/L, respectively, as reported in Table 3. Nevertheless, in the IDW map, there was no high <sup>228</sup>Ra content predicted in the easternmost areas of the Ciudad Rodrigo-Salamanca hydrogeological unit, thus groundwaters of this area must be affected by other local conditions (i.e., groundwater chemical composition and/or <sup>232</sup>Th concentration in the rock mineral-bearing) that led to minimizing the presence of <sup>228</sup>Ra.

According to the IDW method (Figure 4e), and the concentration values from the analysed groundwater intakes, <sup>210</sup>Pb > 200 mBq/L (MRC) was not found in the DB. Significant high concentrations, between 100 and 200 mBq/L, were predicted only in the Ciudad Rodrigo-Salamanca and Los Arenales hydrogeological units, with maximum values exceeding 150 mBq/L at some sampling points, where groundwaters are stored in the Cenozoic and Quaternary medium permeability detritic lithologies. The <sup>210</sup>Pb presence in near-neutral pH groundwaters is mainly related to the alpha recoil of <sup>214</sup>Pb and its solubility enhancement under reducing conditions [25,53]. Moreover, as <sup>210</sup>Pb is a progeny of <sup>226</sup>Ra, its presence may be associated with a high <sup>226</sup>Ra presence in groundwater. The similarities observed between the spatial distribution of <sup>226</sup>Ra and <sup>210</sup>Pb in the Ciudad Rodrigo-Salamanca unit, as illustrated in Figure 4c,e, suggest that the same hydrogeological

settings are the controlling processes involved in the water–rock interaction and are also the same.

#### *3.2. Radioactivity Parameters Used for Drinking Water Monitoring*

As has been discussed in the previous section, the spatial distribution of the uranium isotopes, <sup>238</sup>U and <sup>234</sup>U, had relatively high concentrations, exceeding 100 mBq/L, in wide areas comprising several hydrogeological units of the DB. Both radioisotopes disintegrate by α emission, so, consequently, their occurrence in groundwater contributes to the *a*<sup>α</sup> value. This relation between uranium concentration and *a*<sup>α</sup> spatial distribution gave rise to a similar prediction map, as illustrated in Figures 4a,b and 5a. Prediction areas with values higher than 100 mBq/L occupied a wide extension of the Castilla y León territory, and, in the southern boundaries of the Los Arenales, Segovia and Ciudad Rodrigo-Salamanca hydrogeological units. even exceeded by one magnitude order, with a maximum of 2407, 1847, 1223 mBq/L, respectively (Table 3). As can be observed in Figure 5a, in the southern areas of the Los Arenales region and the Segovia HU, the high concentration seemed to follow the groundwater flow pattern towards the centre of the DB, which suggests transference of α emitters, mainly uranium radionuclides, through the unconfined DB aquifers of high concentrations, from the highly influenced mineralogical detritic aquifers to the Duero River. The shape of the interpolated *a*<sup>α</sup> matched with the predicted areas of high uranium concentration, as illustrated in Figure 4a,b. In the shallow and well-oxygenated aquifers, formed in alluvial and terraces deposits of the current riverine system occupying the central areas of the DB, *a*<sup>α</sup> values were found. Regarding *a*<sup>α</sup> prediction in the central carbonate aquifers, especially significant was the concentration in Paramo de Torozos HU, which was higher than 1100 mBq/L, due to the dissolution of the limestones leading to an increase in the concentration of uranium isotopes in groundwater. This correlation between uranium occurrence in the Spanish carbonate aquifers has been previously demonstrated [10,15].

Regarding the *a*<sup>β</sup> prediction map, the occurrence of higher values than their respective screening level, 1000 mBq/L, was limited to small regions within the DB, as seen in Figure 5b. The long-lived <sup>238</sup>U-daughters, <sup>234</sup>Th and 234mPa, that decay from β-emission, were in secular equilibrium with their parent in the analysed samples. Both β-emitting radionuclides were the main contributors to *a*β, especially in the sandy aquifers of the Esla-Valderaduey region and the Páramo de Torozos HU, with maximum values of 1547 and 1425 mBq/L, respectively. It should be noted that *a*<sup>β</sup> hot spots were also influenced by the contribution of the highest <sup>228</sup>Ra and <sup>210</sup>Pb concentrations, such as those occurring in groundwaters from the southwestern of the Ciudad Rodrigo-Salamanca HU.

According to the *ID* prediction map shown in Figure 4f, the radioactivity content in groundwaters from the southwestern areas of Ciudad Rodrigo-Salamanca HU, as well as the aquifers placed in the southwestern of the Esla-Valderaduey region, predicted *ID* that exceeded their corresponding parametric value. The maximum *ID* values of 0.32 and 0.17 mSv corresponded to both hydrogeological units, respectively, as reported in Table 3. This meant that the consumption of these groundwaters may pose a risk to public health. The suitability of the screening levels as a monitoring strategy is out of the scope of this research, and, due to the low selectivity of screening levels, especially the *a*<sup>α</sup> not really appropriate. The interpolation maps do, however, provide evidence that performing specific radionuclide analysis to properly evaluate the radioactivity quality of the groundwaters for human consumption, is necessary.

**Figure 5.** Spatial distribution of the gross alpha activity (**a**) and gross beta activity (**b**) screening levels for the hydraulically unconfined aquifer system of the Duero Basin. The hydrogeological unit names corresponding to the numbers illustrated are referred to in Figure 2. **Figure 5.** Spatial distribution of the gross alpha activity (**a**) and gross beta activity (**b**) screening levels for the hydraulically unconfined aquifer system of the Duero Basin. The hydrogeological unit names corresponding to the numbers illustrated are referred to in Figure 2.

#### **4. Conclusions**

As part of more extensive research concerning the groundwater radiological characterisation of CyL, this paper addressed the use of a proper interpolation method, IDW, for elaborating radioactivity prediction maps to provide an overview concerning radioactive spatial distribution in the hydraulically interconnected aquifers of the DHC. The major occurrence of <sup>234</sup>U and <sup>238</sup>U activity concentrations were found in the southern area of the DB, mainly promoted by the groundwater flux toward the Duero River, leaching the minerals of the dentritic materials, which form a moderate to highly productive aquifer system, the composition of which is influenced by the igneous and metamorphosed rocks of the CS. Despite the presence of the radionuclides from the natural U- and Th-decay series in groundwater being controlled by several interconnected hydrogeochemical and physical parameters, such as redox, pH, adsorption/desorption rates, water composition, alpha recoil, etc., and water-rock interaction, the geostatistical analysis performed showed that the radioactivity content in DHC groundwaters is related to the hydrogeological context.

It should be underlined that further hydrogeological evaluation is necessary for a better understanding of the behaviour of radionuclides in groundwaters, which was out of the scope of this study. The prediction maps generated a worthwhile qualitative overview concerning the spatial distribution of radionuclides in the studied area. Maps can contribute towards improving radioactivity quality control planning, and in identifying hazardousprone areas, in which the consumption of groundwaters may pose a potential risk to human health. According to the current drinking water framework, and given the radionuclide spatial distribution content in the DB, the use of the prediction maps of *a*<sup>α</sup> and *a*<sup>β</sup> may not be suitable as a radioactivity screening strategy for cost-effective management in groundwater quality control.

**Author Contributions:** Conceptualization, D.B.-A.; methodology D.B.-A., B.Q., A.M.M.-G. and J.C.L.; software, D.B.-A.; validation, D.B.-A., B.Q., A.M.M.-G. and J.C.L.; formal analysis, D.B.-A.; investigation, D.B.-A. and B.Q.; resources, D.B.-A. and B.Q.; data curation, D.B.-A., B.Q. and A.M.M.-G.; writing—original draft preparation, D.B.-A.; writing—review and editing, D.B.-A.; visualization, D.B.-A.; supervision, D.B.-A., A.M.M.-G., B.Q. and J.C.L.; project administration, B.Q.; funding acquisition, B.Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors wish to acknowledge the staff of the Dirección General de Salud Pública of Junta de Castilla y León for the planning, sampling and discussion of the results during the course of the project. This work was promoted and funded by the Dirección General de Salud Pública (Consejería de Sanidad de la Junta de Castilla y León), in the framework of the research project "Radiological characterisation of groundwater intakes intended for human consumption".

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**

