A European map of groundwater pH and calcium

Abstract

Water resources and associated ecosystems are becoming highly endangered due to ongoing global environmental changes. Spatial ecological modelling is a promising toolbox for understanding the past, present and future distribution and diversity patterns in groundwater-dependent ecosystems, such as fens, springs, streams, reed beds or wet grasslands. Still, the lack of detailed water chemistry maps prevents the use of reasonable models to be applied on continental and global scales. Being major determinants of biological composition and diversity of groundwater-dependent ecosystems, groundwater pH and calcium are of utmost importance. Here we developed an up-to-date European map of groundwater pH and Ca, based on 7577 measurements of near-surface groundwater pH and calcium distributed across Europe. In comparison to the existing European groundwater maps, we included several times more sites, especially in the regions rich in spring and fen habitats, and filled the apparent gaps in eastern and southeastern Europe. We used random forest models and regression kriging to create continuous maps of water pH and calcium at the continental scale, which is freely available also as a raster map (Hájek et al., 2020b; 10.5281/zenodo.4139912). Lithology had a higher importance than climate for both pH and calcium. The previously recognised latitudinal and altitudinal gradients were rediscovered with much refined regional patterns, as associated with bedrock variation. For ecological models of distribution and diversity of many terrestrial ecosystems, our new map based on field groundwater measurements is more suitable than maps of soil pH, which mirror not only bedrock chemistry but also vegetation-dependent soil processes.

1 Introduction

The Earth system is currently undergoing unprecedented changes in climate, global biogeochemical cycles and land use, resulting in biodiversity loss (Ceballos et al., 2017; Song et al., 2018; Blowes et al., 2019; Brondizio et al., 2019). Freshwater systems belong to the most endangered habitats (Cantonati et al., 2020a; Tickner et al., 2020), and, among them, groundwater-dependent ecosystems, such as fens and springs, hold primacy (Janssen et al., 2016; Chytrý et al., 2019; Hájek et al., 2020a; Stevens et al., 2020). Species composition and richness of spring systems are generally governed by water pH and calcium concentration (${\mathrm{Ca}}^{\mathrm{2}+}$), which are highly variable at different spatial scales (Malmer, 1986; Rydin et al., 2013; Peterka et al., 2017; Horsáková et al., 2018; Cantonati et al., 2020a, b). Therefore, understanding the spatial patterns in groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ is important not only for general geochemical knowledge and for water resource management but also to the same extent for the conservation of freshwater systems and associated biodiversity.

In Earth and biodiversity sciences, ecological modelling is a widely used tool for understanding the distribution and diversity patterns of ecosystems and habitats and for predicting their future development under global change. Ecological models usually incorporate environmental or historical predictors extracted from thematic maps (Jiménez-Alfaro et al., 2018a; Večeřa et al., 2019; Divíšek et al., 2020), including soil properties for terrestrial ecosystems (Hengl et al., 2017). However, soil parameters such as soil pH contribute negligibly to the models for groundwater-dependent habitats, even for those strongly controlled by pH and ${\mathrm{Ca}}^{\mathrm{2}+}$, such as base-rich fens (Jiménez-Alfaro et al., 2018b). This is due to a poor correlation between groundwater chemistry and pH or ${\mathrm{Ca}}^{\mathrm{2}+}$ in soil, disrupted mainly by mineral leaching or accumulation of organic matter in soil. For this reason, there is a strong need to produce maps for groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration at the European scale that would allow for producing the continental-scale ecological models useful for enforcing conservation strategies in groundwater-dependent habitats. Ideally, such models should include lithology as a dominant factor determining groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration (Hem, 1985; Chapelle, 2003; Tahvanainen; 2004; Stevens et al., 2020).

In spite of important mapping efforts of groundwater (Duscher et al., 2015) and karst aquifers (Chen et al., 2017) at the European and global level, the only available European-scale maps of groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration are those included in the FOREGS (Forum of European Geological Surveys) Geochemical Atlas of Europe (Salminen et al., 2006). These maps are based on 808 stream-water measurements distributed relatively equally across Europe. However, they show a large gap in eastern and southeastern Europe (Romania, Bulgaria, Belarus, the Russian Federation, Ukraine, Moldova, Serbia, Kosovo, Montenegro, Bosnia and Herzegovina, and North Macedonia). In addition, those maps are based on insufficient data density in some areas rich in groundwater-dependent ecosystem types but heterogeneous in terms of lithology (the Alps, the Carpathians, Bohemian Massif, the Cantabrian Mountains and the Pyrenees, and some regions of Fennoscandia). We therefore aimed at substantial improvement of the existing data by creating a database with field data measurements across the entire European continent and at creating a model-based map representing major patterns of groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration at local and continental scales. Our data will allow for a better understanding of the patterns and causes of groundwater conditions in freshwater systems, strongly improving the spatial information suitable for European-scale modelling of biodiversity in groundwater-dependent and related ecosystems.

Table 1

Data sources. The name abbreviations are explained in the team list or Acknowledgements. LANUV: Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen (State Agency for Nature, Environment and Consumer Protection North Rhine-Westphalia).

(I) Vegetation databases
Name	$n$	GIVD code	Custodian
European Mire Vegetation Database (EMVD)	510	GIVD EU-00-022	T.P.
National Vegetation Database of Denmark	373	EU-DK-002	J.E.M.
UK National Vegetation Classification Database	224	GIVD EU-GB-001	J.R.
Balkan Vegetation Database	27	GIVD EU-00-019	K.V.
The vegetation database of Mecklenburg-Vorpommern	19	GIVD EU-DE-001	F.J.
Basque Country	7	EU-00-011	I.B.
(II) Unpublished data sets
Regions	$n$	Co-authors of the data set
Central and eastern Europe	1405	Z.P., M.Há., P.H., T.P., D.D., L.S., M.L., J.N., P.P., P.S., A.Š., Y.S., C.B.-N.
Spain	645	A.P-H., E.P-I., B.J-A.
Bulgaria	428	P.H., M.Há., M.Hor.
Fennoscandia	392	M.Há., T.P., D.D., M.Hor., V.H., P.H., T.K., T.T., J.Kap., D.-I.Ø.
Apennines	285	M.T., M.Can., M.Car., S.S., A.P., L.B., R.G.
Europe (cross-taxon research)	281	M.Há., P.H., D.D., M.Hor., V.H.
Balkans except Bulgaria	134	A.D., E.M., J.Kam., P.L., T.P., M.Há., P.H.
(III) Public data sets
Area and agency	$n$	Provided via
LANUV (North Rhine-Westphalia, Germany)	463	Dr. Dirk Hinterlang, Dr. Sabine Bergmann
Ministry of Environment and Water (Bulgaria)	442	Mrs. Rossitza Gorova
Ministry of the Environment and Spatial Planning (Slovenia)	176	http://www.arso.gov.si/ (last access: 26 February 2019)
(IV) Geochemical Atlas of Europe
Area	$n$	Reference
Europe	808	Salminen et al. (2006)
(V) Other literature data (gap-oriented excerption)
Region and context	$n$	Reference
West-central European springs	340	Hinterlang (1992)
Lithuanian springs	194	Kadūnas et al. (2017)
Northeastern England	111	Eades et al. (2018)
Northern Germany	58	Tanneberger et al. (2011)
Central Bosnia	50	Savić et al. (2017)
Western Bohemian mineral springs (Czech Republic)	28	Laburdová and Hájek (2014)
Scotland	24	Gorham (1957)
Eastern Bosnia	20	Kamberović et al. (2019)
Switzerland (mires)	14	Lamentowicz et al. (2010)
British and French travertines	13	Pentecost and Zhaohui (2002)
Northwestern Poland (mires)	8	Lamentowicz and Mitchell (2005)
Western Balkans	6	Ridl et al. (2018)
Kosovo	5	Kelmendi et al. (2018)
Southeastern Croatia	4	Terzić et al. (2014)
Kosovo (Rugova)	3	Lajçi et al. (2017)
Serbia	3	Ćirić et al. (2018)
Northeastern Croatia	1	Špoljar et al. (2011)
Southeastern Croatia (Krčić)	1	Kolda et al. (2019)

2 Methods

2.1 Data collection

We assembled the data set of pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ (or electrical conductivity in $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$ at 20 ${}^{\circ }\mathrm{C}$; hereinafter abbreviated as EC) measurements in groundwater, covering the whole of Europe, with a greater density in the regions rich in endangered groundwater-dependent ecosystems such as springs and fens. We excluded most of Ukraine and the European part of the Russian Federation because of large data gaps in these areas. We considered all types of shallow groundwater systems, especially spring, spring fen and stream water. The core of our data set is formed by unpublished pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ or EC data sets of the co-authors (3618 sites); some of them were processed in ecological papers without having presented original pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ data (Hájková et al., 2006, 2008; Hájek et al., 2008; Sekulová et al., 2013; Plesková et al., 2016; Horsáková et al., 2018; Šímová et al., 2019). The second most important source were vegetation databases registered in GIVD (Global Index of Vegetation‐Plot Databases; Dengler et al., 2011; Table 1) and EVA (European Vegetation Archive; Chytrý et al., 2016), from where 1160 measurements from freshwater habitats were obtained. Both unpublished data and data from vegetation databases were filtered using original information or metadata of the sources in a way that only data from spring-fed fens and springs were considered. The data from ombrotrophic bogs and clearly topogenic fens (mainly terrestrialised lakes) were omitted because their water chemistry is governed by the decomposition of organic matter, atmospheric humidity and deposition, algal photosynthesis (Kann and Smith, 1999), and biotic processes such as cation exchange capacity of mosses (Clymo, 1963; Soudzilovskaia et al., 2010; Vicherová et al., 2015) rather than by bedrock chemistry. We also obtained data from public data sets stored in national environmental and nature conservation agencies of Germany, Slovenia and Bulgaria (1081 sites; see Table 1); data from the FOREGS Geochemical Atlas of Europe (Salminen et al., 2006; 808 sites); and literature data based on our gap-oriented excerption (883 sites; Table 1); most data came from Hinterlang (1992), Tanneberger et al. (2011), Eades et al. (2018), Kadūnas et al. (2017) and Savić et al. (2017).

In total, we collected 7577 samples (Table 1). Some of these samples are repeated measurements conducted in the same site, especially in public data sets, while other samples (from vegetation databases or literature data) share the same coordinates and site name or code, suggesting repeated measurements as well. We therefore averaged repeated measurements from the same sampling sites. We further deleted samples whose coordinates were obviously erroneous, such as those in oceans. These steps reduced the number of samples to 6561, out of which 6459 samples contained information on the water pH, and 5927 samples contained information about EC of water or ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration. Out of these 5927 samples, 2988 had directly measured both ${\mathrm{Ca}}^{\mathrm{2}+}$ and EC ($\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$ at 20 ${}^{\circ }\mathrm{C}$), and for the remaining 2939 samples we estimated the ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration by EC of water.

Figure 1

The final regression model to impute ${\mathrm{Ca}}^{\mathrm{2}+}$ values based on electrical conductivity (EC; $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$; $n$ $=$ 2062).

[Figure omitted. See PDF]

2.2

Imputation of missing ${\mathrm{Ca}}^{\mathrm{2}+}$ values by EC of water

For imputation of ${\mathrm{Ca}}^{\mathrm{2}+}$ values based on EC, we first aimed at constructing a simple imputation equation based on the well-known correlation between EC and the ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration in springs and fens (Hem, 1985; Sjörs and Gunnarsson, 2002; Plesková et al., 2016). In our data set of 2988 samples, as well as in its regional subsets, this relationship was strongly governed by EC values above ca. 1000 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$, although they formed only a small part of the data set (4.7 % of the data set; 139 samples). In the EC range of 1000–10 000 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$ (an outlier of 17 000 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$ was omitted), the correlation between water EC and the ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration was not statistically significant ($r$ $=$ 0.15, $p$ $=$ 0.07). The problem of high EC values governing the regression model was the most apparent in the public data sets. In the data set of the Bulgarian Ministry of Environment and Water, weak correlation between EC and ${\mathrm{Ca}}^{\mathrm{2}+}$ persisted even when EC values above 1000 were omitted (Fig. S1 in the Supplement). This database further contains many samples which are not near-surface samples that were measured in other datasets. We therefore finally decided (1) not to include the database of the Bulgarian Ministry of Environment and Water into the imputation model and (2) limit the gradient of EC to 1000 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$. We further omitted a few samples from ophiolite (Kamberović et al., 2019) where high EC occurred despite low ${\mathrm{Ca}}^{\mathrm{2}+}$. The resulting data set of 2319 samples nevertheless still showed some samples with a suspiciously high or low ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration relative to EC (Fig. S2 in the Supplement), suggesting either the effect of ions other than ${\mathrm{Ca}}^{\mathrm{2}+}$ or an inconsistent analytical methodology. Because our aim was to create the most accurate imputation model rather than test the relationship, we removed these outliers. Therefore, we calculated the $\mathrm{EC}:\mathrm{Ca}$ and $\mathrm{Ca}:\mathrm{EC}$ ratios and removed outliers, i.e. all points larger than 1.5 times the interquartile range. The final imputation model was hence based on 2062 sites. We performed a null-intercept linear regression (Fig. 1) with ${\mathrm{Ca}}^{\mathrm{2}+}$ as a dependent variable ($y$) and EC as a predictor ($x$); the resulting equation of $y$ $=$ 0.153$x$ was obtained ($R^{\mathrm{2}}$ $=$ 0.84). Such a relationship between ${\mathrm{Ca}}^{\mathrm{2}+}$ and EC is similar to that found in the abovementioned studies (Hem, 1985; Sjörs and Gunnarsson, 2002; Plesková et al., 2016). Based on this equation, we imputed ${\mathrm{Ca}}^{\mathrm{2}+}$ concentrations to all samples where only EC was measured. The imputed ${\mathrm{Ca}}^{\mathrm{2}+}$ values show a somewhat narrower range (Fig. S3 in the Supplement) than originally measured values. Both subsets show a minimum ${\mathrm{Ca}}^{\mathrm{2}+}$ value below 1 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}$, but imputed data show a lower non-outlier maximum (125.5 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}$) than measured data (197.2 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}$). The absolute maximum value was also lower for the subset with imputed values. Imputation of ${\mathrm{Ca}}^{\mathrm{2}+}$ values to all samples, including high EC ones ($>$ 1000 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$), hence did not skew the imputed data to higher values.

2.3 Geographical modelling and selection of the predictors

We used our dataset with measured groundwater pH and either measured or imputed ${\mathrm{Ca}}^{\mathrm{2}+}$ concentrations to model expected values across non-sampled areas. Our aim was to produce continuous maps for groundwater-dependent pH (GW-pH) using 6459 samples (pH min $=$ 2.20; max $=$ 11.32; mean $=$ 6.69) and groundwater-dependent ${\mathrm{Ca}}^{\mathrm{2}+}$ (GW-Ca) using 5927 samples (min $=$ 0.15; max $=$ 3567.41; mean $=$ 48.73 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}$). ${\mathrm{Ca}}^{\mathrm{2}+}$ values were ln-transformed. All field samples had geographic coordinates assigned from GPS or georeferenced with an accuracy between ca. 10 $\mathrm{m}$ (precise field measurements) and 500 $\mathrm{m}$ (from georeferenced sites in databases). We kept the pH outliers: 10 values below 3.5 and 9 values above 8.8. Even if these values may be suspicious, they largely come from published sources (FOREGS Geochemical Atlas of Europe and the British vegetation database). Apart from measurement error, they may be explained by the influence of mineral waters from deep hydrological circulations (e.g. in a spring in the Apennines, a very high pH value of 11.2 was due to enrichment with sodium and chloride associated with a low-temperature reaction between meteoric water and ultramafic rocks; Boschetti and Toscani, 2008; Boschetti et al., 2013; Segadelli et al., 2017; Cantonati et al., 2020c). These values form only a minor part of the data set and have a negligible effect on the results.

For each site, we obtained environmental predictors from thematic GIS (geographic information system) maps (see below). We focused on the predictors that may causally affect the groundwater pH and calcium concentration. Aquifer chemistry is of prime importance (Hem, 1985; Fairchild et al., 1994; Frei et al., 2000; Chapelle, 2003; Tahvanainen, 2004; Stevens et al., 2020), but no such thematic map exists, at least not at the scale needed for computing our spatial predictions. We therefore included the lithological groups from the Hydrogeological Map of Europe (Duscher et al., 2015), together with soil pH maps (see below), for which we anticipated a certain correlation with bedrock chemistry. Apart from aquifer chemistry, residence time may also affect groundwater chemistry by impacting dissolution rates. Precipitation amount and frequency affect not only flow paths activity and redistribution of groundwater but also its residence time in the aquifer, impacting carbonate dissolution and precipitation rates (Hem, 1985; Crossman et al., 2011; Lewandowski et al., 2015; Vystavna et al., 2020). Groundwater with a short transit time (1–3 years) or “young water” (Soulsby et al., 2015) can be particularly sensitive to changes in precipitation amount and frequency. We therefore also considered climatic parameters associated with precipitation in the models (see below). Although there are some other potential predictors of minor importance that may affect groundwater chemistry (Hem, 1985; Stevens et al., 2020), no corresponding thematic map is available to be included into our models. For some sites, the selected predictors were missing in the maps (e.g. sites at the far north in the Arctic zone or close to the sea or water bodies), and these sites were therefore not included in the final models. We finally collected topographic data to test the potential effect of elevation and slope as indirect factors potentially influencing groundwater chemistry.

Figure 2

Spatial distribution of the three groups of calibration data collected for modelling groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ in European fens (original and literature data from springs and fens; data from streams from the FOREGS Geochemical Atlas of Europe; other data). Other data include public data from national groundwater monitoring in Bulgaria and Slovenia. For separate maps of pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ see Fig. S5 in the Supplement.

[Figure omitted. See PDF]

2.4 Numerical analyses

Numerical analyses were done in R version 3.6.3 (R Core Team, 2020), with the support of ArcGIS 10.2 (Esri, Redlands, CA, USA) for geoprocessing and map production. We first conducted exploratory analyses to test the prediction ability of GIS layers related to soil bedrock, climate and topography on the variation of both Ca and pH. We focused on layers with a complete coverage of Europe, with an eastern border of the Black Sea in Turkey to the White Sea in the Russian Federation, thus including the regions with a relatively good cover of field measurements (Fig. 2). We performed linear models for individual variables to select those providing significant relationships and $>$ 1 % of explained variance. A variable for soil pH (measured in water solution) at 15 $\mathrm{cm}$ depth for a 250 $\mathrm{m}$ grid resolution provided by the SoilGrids project (https://www.soilgrids.org, last access: June 2020) had the highest explanatory power for GW-pH ($R^{\mathrm{2}}$ $=$ 0.22) and GW-Ca ($R^{\mathrm{2}}$ $=$ 0.16). The same results were obtained when using the same variable for 5 or 10 $\mathrm{cm}$ depth. We also tested soil estimates from Ballabio et al. (2019), but they provided weaker relationships for both GW-pH ($R^{\mathrm{2}}$ $=$ 0.14 using soil pH as a predictor) and GW-Ca ($R^{\mathrm{2}}$ $=$ 0.01 using soil pH; $R^{\mathrm{2}}$ $=$ 0.01 using soil ${\mathrm{CaCO}}_{\mathrm{3}}$). To account for lithology, we used the lithological groups (litho3 level, identifying 12 major classes; see Fig. 4) included in the polygon layer of the Hydrogeological Map of Europe (Duscher et al., 2015) as a categorical variable. We also selected annual precipitation (coded as Bio12) as provided in CHELSA (Climatologies at high resolution for the earth's land surface areas; Karger et al., 2017) to account for precipitation gradients which are expected to influence groundwater regimes. Other CHELSA variables related to precipitation were highly correlated with annual precipitation (Pearson $r$ $>$ 0.75) and omitted. Slope and elevation showed negligible effects on both GW-pH and GW-Ca (linear regression, adjR2 $<$ 0.05, $p$ $<$ 0.001). Since preliminary models showed no differences when these variables were included, they were discharged.

The variables of lithology and soil pH were aggregated to the same grid extent of CHELSA at 1 $\mathrm{km}$ resolution, as this is the most appropriate scale to balance the original scales of both layers. This grid extent is also the most suitable spatial scale to be used in the context of further ecological modelling, which is in many cases combined with climatic data from e.g. CHELSA or WorldClim (http://www.worldclim.org, last access: June 2020) for making temporal climatic projections. The lithological map (originally at 1 : 1 500 000 scale, which corresponds to a raster resolution of ca. 1 $\mathrm{km}$) was converted to a grid resolution using the dominant unit. Soil pH was converted from the original 250 $\mathrm{m}$ to 1 $\mathrm{km}$ grid resolution using a bilinear interpolation to create a smooth surface based on the weighted average of the four nearest cells.

Because even the combination of the selected variables might not lead to a precise fine-scale indication of aquifer chemistry, we further employed a kriging approach to data analysis, assuming an effect of spatial correlation to estimate the values close to the original samples. The final spatial predictions were therefore based on regression kriging (RK), a technique that combines a regression model based on explanatory variables with the interpolation of model residuals with ordinary kriging (Hengl et al., 2007; Meng et al., 2013). RK is especially appropriate for modelling soil attributes at medium and large scales, combining the spatial autocorrelation of soil variables with the explanatory power of auxiliary variables (Keskin and Grunwald, 2018). We implemented RK with the “GSIF” R package (Hengl, 2020). As the regression component, we computed random forests, since a preliminary analysis with our data showed better performance than linear models, generalised linear models or generalised additive models. Random forests are ensemble learning methods based on decision trees and an internal correction of overfitting, which provide high interpretability and good performance when compared with other algorithms used in soil spatial modelling (Wiesmeier et al., 2011). Another advantage of random forests is that they have no requirements for considering the probability distribution of soil variables, fitting complex non-linear relationships for spatial extrapolation (Hengl et al., 2015) that ultimately improve the spatial predictions. We fitted the random forest model and the residual variogram for groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ separately using the function fit.gstatModel() in the “GSIF” package. Effect plots for the predictors were created for the same models using the partial() function in the “pdp” package (Greenwell, 2017). Spatial predictions were then computed with the predict() function using the model object generated previously and a 5-fold cross-validation. Model evaluation was based on the calculation of the mean error (ME) and the root mean squared error (RMSE) as the differences between predicted and observed values (Keskin and Grunwald, 2018; Pham et al., 2019). We compared the relationships between the models produced for both groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ by using a random sampling of 5000 points to extract cell values and computing a Pearson correlation. To assess regional differences, we correlated values grouped in 25 neighbouring cells of each single cell using the rasterCorrelation() function in the “spatialEco” R package (Evans, 2020).

Figure 3

Variable importance of random forest models computed for groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$. MSE: mean standard error.

[Figure omitted. See PDF]

Figure 4

Partial dependence plots showing the effects of the variables used in the random forest models computed for model groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ in Europe.

[Figure omitted. See PDF]

3 Results

In measured data, ranges and medians of pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration were similar across Europe (Fig. S4 in the Supplement), with the lowest pH values found in the Atlantic and Iberian regions and the highest pH values found in southern Europe except the Iberian Peninsula. The lowest ${\mathrm{Ca}}^{\mathrm{2}+}$ values were found in boreal Europe, while the highest were found in central and southern Europe. The random forest models computed with the lithology, soil pH and precipitation explained 40 % and 55 % of the variance for GW-pH and GW-Ca, respectively. Lithology was the variable with the highest importance in both models (Fig. 3), although its effect was higher in the model computed for ${\mathrm{Ca}}^{\mathrm{2}+}$ than for pH. These effects were mainly associated with the lithological units reflecting calcareous bedrock, followed by categories with coarse and fine sediments such as flysch (Fig. 4). Soil pH had a higher relative importance in GW-pH than GW-Ca, although in both cases the variable had a similar positive effect. Finally, annual precipitation had the lowest contributions in the two models (Fig. 3), with both ${\mathrm{Ca}}^{\mathrm{2}+}$ and pH dropping suddenly after the threshold of annual precipitation of ca. 1800 $\mathrm{mm}$ (Fig. 4), although the highest pH values occur under the lowest precipitation and tend to decrease towards high-precipitation areas.

When adding the kriging component, model predictions reached 65 % and 74 % of explained variance for GW-pH and GW-Ca, respectively. The mean values of standard errors (SEs; 0.0006 for pH; $-$0.0009 for Ca) and root mean squared errors (RMSEs; 0.588 for pH; 0.690 for Ca) were higher in the models for pH, but in both cases they showed low values that suggest accurate predictions, in agreement with their total explained variance.

Figure 5

Model predictions based on regression kriging. Note that the ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration is on an ln scale. The map was created entirely by the authors with state borders from ArcGIS® version 10.2 (Esri, Redlands, CA, USA).

[Figure omitted. See PDF]

Figure 6

Spatial correlation between the models computed for groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$. Values show the Pearson correlation coefficient computed over every single cell by using a sampling of 25 neighbouring cells. The map was created entirely by the authors with state borders from ArcGIS® version 10.2 (Esri, Redlands, CA, USA).

[Figure omitted. See PDF]

Model predictions for GW-pH reflected the lowest values in Scandinavia, Iceland, the northern UK, and some regions of central and eastern Europe (Fig. 5). The highest values were predicted in eastern Iberia and many regions of central and eastern Europe, although a big part of the study area was dominated by neutral pH values (6 to 7). The spatial patterns for GW-Ca (Fig. 5) were rather similar to pH. The overall correlation between the two models was 0.83 (Pearson $r$, $p$ $<$ 0.001), but they showed differences in some regions. This was supported by the spatial correlation computed for each cell (Fig. 6), reflecting different magnitudes of correlation across the study area, especially in the eastern Iberian Peninsula and southeastern Europe.

4 Discussion

4.1

Spatial patterns in groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration in Europe

As expected, the values of water pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration are largely shaped by lithology in groundwater-dependent habitats across Europe. Indeed, it has been recognised by regional studies that the distribution of major spring and fen habitats, of which the species composition largely depends on pH and ${\mathrm{Ca}}^{\mathrm{2}+}$, is well determined by bedrock type (Hájek et al., 2002; Tahvanainen, 2004; Hinterlang, 2017; Peterka et al., 2017; Cantonati et al., 2020c). Since the European-scale geological map we used here is not precise enough to capture differences in bedrock chemistry within the major lithological units that are defined largely by geological stratification, the contribution of soil pH in the model probably also reflected lithological variation, as soil pH generally correlates with regional bedrock chemistry (Chadwick and Chorover, 2001). On the other hand, soil pH is also affected by climate-dependent pedogenesis, which incorporates a climate-zonal geographical component in this effect (Duchaufour, 2012; Maxbauer et al., 2017).

Precipitation is another determinant of groundwater chemistry in our study. High annual precipitation above ca. 1800 $\mathrm{mm}$ obviously reduced an interaction time of groundwater with ${\mathrm{Ca}}^{\mathrm{2}+}$ and carbonates deposited in rocks (Fairchild et al., 1994; Segadelli et al., 2017; Cantonati et al., 2020b), resulting in lower ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration in groundwater. This effect is more pronounced in snowy regions, where seasonal snowmelt modulates the recharge patterns of groundwater. The duration of the snowmelt period can impact the occurrence and dynamic of preferential flow and prolong or reduce the interaction of the seepage with soil and bedrock materials (Mohammed et al., 2019). We therefore suggest that fast hydrological pathways and a short transit time driven by snowmelt and precipitation can explain the lowest ${\mathrm{Ca}}^{\mathrm{2}+}$ values in hyper-oceanic cold regions of southwestern Norway or western Scotland. It may further explain lower pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ values on windward slopes of high mountains, even if bedrock is moderately calcium-rich. Nevertheless, understanding the complementary effects of precipitation and slope will need to account for more accurate models based on GPS data, better precipitation data and high-resolution ($<$ 250 $\mathrm{m}$) topographic predictors.

The resulting pattern at the European scale is governed by the strong latitudinal and altitudinal gradients, i.e. decreasing pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ northwards and regionally also towards mountain regions. This pattern largely follows bedrock chemistry, with crystalline rocks prevailing and most carbonate rocks being eroded by glaciers, in high latitudes and altitudes. The excess of precipitation over evaporation and theoretically also slower weathering rates in colder regions (White and Blum, 1995) contribute as well. Although this pattern is well known (Økland et al., 2001; Hájek et al., 2006; Hinterlang, 2017; Peterka et al., 2017) and has been captured also by the FOREGS Geochemical Atlas of Europe (Salminen et al., 2006), our improved model provides much finer regional patterns. In southern Europe, low pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ values were modelled in the Pyrenees, the Balkans, southwestern Corsica and Calabria, i.e. the regions where boreal or endemic types of fen communities occur as relicts (Chytrý et al., 2020). The Alps, the Apennines, the Carpathians and the Baltic region show a fine-scaled mosaic of alkaline (calcium-rich) and acidic (calcium-poor) groundwater that contributes to the high diversity and conservation value of groundwater-dependent ecosystems, such as fens (Cantonati et al., 2009, 2011; Gerdol et al., 2011; Joosten et al., 2017; Horsáková et al., 2018). The most apparent “acidic island” in central Europe is located in the southwestern part of the Bohemian Massif (Czech Republic and Germany), where acidic types of springs and fens are quite frequent, and some studies further document anthropogenic acidification on siliceous bedrock in the 1970s and 1980s, which has re-emerged recently because of extreme climatic events (Kapfer et al., 2012; Schweiger et al., 2015). It is, however, possible that particularly this acidic island is picked out mainly because of the high amount of available data.

Clearly, most of Fennoscandia is markedly acidic and calcium-poor mainly due to glacial history. Yet, the model identified small alkaline- and calcium-enriched islands in northeastern and central Sweden and northwestern Norway which are associated with rare types of calcareous fen and spring communities (Dierssen, 1982; Vorren et al., 1999; Udd et al., 2015; Miller et al., 2021). More localised pockets of calcareous habitats are however known from most parts of Fennoscandia that are not recognised with the grain of our European-wide analysis. With our results, the future modelling of diversity and distribution of individual habitat types of groundwater-dependent wetlands will be more reliable. Regionally rare habitat conditions will be recognised better, and the disentangling of climate and pH effects will be more feasible.

4.2 Data gaps and further improvements

Although being based on the hitherto most comprehensive field data set currently available, the presented map cannot be considered definitive. Surely there are many pH and EC or ${\mathrm{Ca}}^{\mathrm{2}+}$ measurements conducted across Europe that we could not include into the data set because they are hardly accessible. Except for the Russian Federation and Moldova, the largest gaps still occur in the southern parts of the Pannonian Plain (southeastern Hungary, northern Serbia and western Romania), southeastern Belarus and eastern Ukraine. We have available some data from the latter region (Vystavna et al., 2015; Table S2 in the Supplement), but a large gap in the rest of the data set prevented reliable geospatial modelling. These data might be used in future updates of the map once the gap in central Ukraine is filled. The lack of data in the Pannonian Plain has led to a poor correlation between predicted pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ values (Fig. 4). Such a poor correlation and sometimes low density of data apply also for some other lowland regions, such as the Danube Plain in southern Romania; Po Valley in Italy; and valleys around the Douro, Ebro and Tagus rivers in Spain. Apart from eutrophication, this pattern may be caused by the imbalanced distribution of groundwater-dependent habitat types in our data set. Unlike mountain regions, the data for these lowlands were largely taken from the FOREGS Geochemical Atlas of Europe (Salminen et al., 2006) and national groundwater databases, i.e. largely from stream water. Considering the major purpose of our map (creating new data sources for ecological modelling of fens and springs) these regions are less crucial for biodiversity modelling because they have largely been transformed to arable land or they are too dry. On the other hand, caution is needed when interpreting the maps in an ecological sense. The extremely high pH ($>$ 8) and ${\mathrm{Ca}}^{\mathrm{2}+}$ ($\ln \left[{\mathrm{Ca}}^{\mathrm{2}+}\right]$ $>$ 4; i.e. ${\mathrm{Ca}}^{\mathrm{2}+}$ $>$ 55 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}$) values that occur in lowlands visually govern the map, but for ecological differentiation of groundwater-dependent habitats in Europe the differences within the middle part of the gradient, i.e. between pH 5.5 and 7.0, are much more important (Malmer, 1986; Wheeler and Proctor, 2000; Hájek et al., 2006; Rydin et al., 2013).

Our data set is expected to be amended in the future, as more studies will be published and more data will be available, so further versions will be accessible in the open repository. New data will help to improve predictions for those regions with a relatively lower sampling effort and also those with lithologically heterogeneous landscapes. Future updates of the model may also focus at finer spatial resolutions (e.g. 100 to 250 $\mathrm{m}$), but this will require an increase of the spatial accuracy of the calibration data and the predictor variables, as in some central European areas (Le et al., 2019; Chuman et al., 2019). Although we tested several variables with a potential predictive effect on groundwater pH and Ca, many of them had lower explanatory power (e.g. Ballabio et al., 2019), or they were redundant with the soil pH layer we used (i.e. previous versions of SoilGrids using the same data sources). The low predictive value of other predictors with potential predictive value, like slope, elevation or another precipitation variable, is probably related to their broad spatial scale, with values averaged at 1 $\mathrm{km}$ grid resolution having little impact to discriminate groundwater variation at the landscape level (Jiménez-Alfaro et al., 2018b). This contrasts with the predictive value of the lithological layer, which is however based on a similar spatial (but originally vectorial) resolution. However, lithological bedrock is generally more homogeneous at the landscape scale, with the only exception of certain geologically complex regions. These findings suggest that future improvements of our models will depend on the quality of new lithological (or related soil chemistry) variables with a direct effect on groundwater pH and Ca. We also note that the lithological map we used here is simplified to large units, while many regional and national geological maps are being produced at finer resolutions. The combination of such new predictors with accurate calibration data at the continental level and at fine-scale resolution seems the most likely opportunity to produce significantly better models, since the influence of spatially dependent (kriging) effects is rather limited by the distribution of sample points.

5 Data availability

The dataset of georeferenced pH and EC measurements and the resulting maps in GIS-compatible format (shapefile) are accessible at https://www.zenodo.org/ (last access: March 2021); 10.5281/zenodo.4139912 (Hájek et al., 2020b).

6 Code availability

No original R code was used; the used codes are cited.

7 Conclusions

Here, we provide the first European map of groundwater pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ content. We collected field measurements, distributed as evenly as possible, of water pH and ${\mathrm{Ca}}^{\mathrm{2}+}$ or EC from European groundwater-dependent habitats, with a relatively higher sampling effort in regions rich in endangered groundwater-dependent ecosystems (springs and fens). Despite the general high accuracy of our models, we note that prediction uncertainties may affect the reliability of models computed with both random forests and kriging (Hengl et al., 2018; Szatmári and Pásztor, 2019). Another source of prediction uncertainty is related to the quality of the original chemical measurements and the georeferentiation of their geographic position. Moreover, the predictor variables rely on spatial models (soil pH and precipitation) or broad geographic maps (lithology) which are based on their own uncertainties and assumptions. Future improvements of groundwater pH and Ca estimates should therefore consider a more accurate set of response variables and fine-scale predictors, preferably including lithology and soil pH (i.e. variables surrogating bedrock chemistry) and precipitation sum (i.e. residence time of groundwater). Despite potential uncertainties and data gaps, this study uses an unprecedented combination of data to provide freely accessible and realistic maps that can be used in any kind of spatial modelling, showing a better resolution and fewer gaps than previously published maps. The character of our input data, which are also freely accessible, predetermines our map to be used in ecological modelling to address the distribution and diversity of groundwater-dependent ecosystems and associated species. We even believe that our maps could be also suitable for ecological modelling of habitats other than those which are groundwater-dependent. It may mirror the bedrock chemistry better than the map of soil pH because soil pH is a resultant of pedogenetic processes, which are tightly associated with the character of the vegetation cover itself (Miles, 1985; Duchaufour, 2012). We conclude that our European maps of near-surface groundwater pH and EC provides the best solution currently available for modelling the biodiversity of groundwater-dependent ecosystems, especially at the continental or supra-regional scale.

The supplement related to this article is available online at: https://doi.org/10.5194/essd-13-1089-2021-supplement.

Author contributions

MH and BJA contributed equally to the paper. They conceived the research, collected data and led the writing process. BJA designed and executed models for both the random forest and regression kriging. OH prepared the input data and final map outputs. All authors provided unpublished data and commented on the paper.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

This work was supported by the Czech Science Foundation (grant no. 19-01775S provided support for Borja Jiménez-Alfaro, Petra Hájková, Veronika Horsáková and Michal Horsák; grant no. GX19-28491X to the Centre for European Vegetation Syntheses provided support for Michal Hájek, Tomáš Peterka and Ondřej Hájek). Jutta Kapfer was supported by the Fram Center (grant no. A36214). We thank Sabine Bergmann and Dirk Hinterlang (State Agency for Nature, Environment and Consumer Protection North Rhine-Westphalia, Germany) and Rossitza Gorova (Ministry of Environment and Water of Bulgaria) for providing us with the public water chemistry data, Andraž Čarni for alerting us to the open data on Slovenian groundwater, and Valerijus Rašomavičius for providing us with the data from Kadūnas et al. (2017). We thank Tatyana Ivchenko for providing us with the data for the Ural Mountains, which were not included in the paper. Paweł Pawlikowski (P.P.), Lucia Sekulová (L.S.), Jana Navrátilová (J.N.) and Dag-Inge Øien (D.-I.Ø.) kindly agreed for us to use their unpublished pH and EC data. We further thank Ilona Knollová, John Rodwell (J.R.), Kiril Vasilev (K.V.) and Idoia Biurrun (I.B.) for providing pH data from vegetation databases via the European Vegetation Archive.

Financial support

This research has been supported by the Czech Science Foundation (grant nos. 19-01775S and GX19-28491X) and the Fram Center (grant no. A36214).

Review statement

This paper was edited by Jens Klump and reviewed by two anonymous referees.