A European map of groundwater pH and calcium

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; https://doi.org/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 (Ca²⁺), 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 Ca²⁺ 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 Ca²⁺, such as base-rich fens (Jiménez-Alfaro et al., 2018b). This is due to a poor correlation between groundwater chemistry and pH or Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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).

2 Methods

2.1 Data collection

We assembled the data set of pH and Ca²⁺ (or electrical conductivity in µS cm⁻¹ at 20 ^∘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 Ca²⁺ or EC data sets of the co-authors (3618 sites); some of them were processed in ecological papers without having presented original pH and Ca²⁺ 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 Ca²⁺ concentration. Out of these 5927 samples, 2988 had directly measured both Ca²⁺ and EC (µS cm⁻¹ at 20 ^∘C), and for the remaining 2939 samples we estimated the Ca²⁺ concentration by EC of water.

Figure 1 The final regression model to impute Ca²⁺ values based on electrical conductivity (EC; µS cm⁻¹; n = 2062).

2.2 Imputation of missing Ca²⁺ values by EC of water

For imputation of Ca²⁺ values based on EC, we first aimed at constructing a simple imputation equation based on the well-known correlation between EC and the Ca²⁺ 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 µS cm⁻¹, 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 µS cm⁻¹ (an outlier of 17 000 µS cm⁻¹ was omitted), the correlation between water EC and the Ca²⁺ 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 Ca²⁺ 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 µS cm⁻¹. We further omitted a few samples from ophiolite (Kamberović et al., 2019) where high EC occurred despite low Ca²⁺. The resulting data set of 2319 samples nevertheless still showed some samples with a suspiciously high or low Ca²⁺ concentration relative to EC (Fig. S2 in the Supplement), suggesting either the effect of ions other than Ca²⁺ 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 EC:Ca and Ca: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 Ca²⁺ as a dependent variable (y) and EC as a predictor (x); the resulting equation of y = 0.153x was obtained (R² = 0.84). Such a relationship between Ca²⁺ 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 Ca²⁺ concentrations to all samples where only EC was measured. The imputed Ca²⁺ values show a somewhat narrower range (Fig. S3 in the Supplement) than originally measured values. Both subsets show a minimum Ca²⁺ value below 1 mg L⁻¹, but imputed data show a lower non-outlier maximum (125.5 mg L⁻¹) than measured data (197.2 mg L⁻¹). The absolute maximum value was also lower for the subset with imputed values. Imputation of Ca²⁺ values to all samples, including high EC ones (> 1000 µS cm⁻¹), 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 Ca²⁺ 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 Ca²⁺ (GW-Ca) using 5927 samples (min = 0.15; max = 3567.41; mean = 48.73 mg L⁻¹). Ca²⁺ values were ln-transformed. All field samples had geographic coordinates assigned from GPS or georeferenced with an accuracy between ca. 10 m (precise field measurements) and 500 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 Ca²⁺ 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 Ca²⁺ see Fig. S5 in the Supplement.

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 cm depth for a 250 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² = 0.22) and GW-Ca (R² = 0.16). The same results were obtained when using the same variable for 5 or 10 cm depth. We also tested soil estimates from Ballabio et al. (2019), but they provided weaker relationships for both GW-pH (R² = 0.14 using soil pH as a predictor) and GW-Ca (R² = 0.01 using soil pH; R² = 0.01 using soil CaCO₃). 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 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 km) was converted to a grid resolution using the dominant unit. Soil pH was converted from the original 250 m to 1 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺. MSE: mean standard error.

Figure 4 Partial dependence plots showing the effects of the variables used in the random forest models computed for model groundwater pH and Ca²⁺ in Europe.

3 Results

In measured data, ranges and medians of pH and Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ and pH dropping suddenly after the threshold of annual precipitation of ca. 1800 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 Ca²⁺ 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 6 Spatial correlation between the models computed for groundwater pH and Ca²⁺. 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).

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 Ca²⁺ concentration in Europe

As expected, the values of water pH and Ca²⁺ 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 Ca²⁺, 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 mm obviously reduced an interaction time of groundwater with Ca²⁺ and carbonates deposited in rocks (Fairchild et al., 1994; Segadelli et al., 2017; Cantonati et al., 2020b), resulting in lower Ca²⁺ 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 Ca²⁺ values in hyper-oceanic cold regions of southwestern Norway or western Scotland. It may further explain lower pH and Ca²⁺ 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 m) topographic predictors.

The resulting pattern at the European scale is governed by the strong latitudinal and altitudinal gradients, i.e. decreasing pH and Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ (ln [Ca²⁺] > 4; i.e. Ca²⁺ > 55 mg L⁻¹) 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 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 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); https://doi.org/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 Ca²⁺ content. We collected field measurements, distributed as evenly as possible, of water pH and Ca²⁺ 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.