Integrating field work and large-scale modeling to improve assessment of karst water resources

Comprehensive management of karst water resources requires sufficient understanding of their dynamics and karst-specific modeling tools. However, the limited availability of observations of karstic groundwater dynamics has been prohibiting the assessment of karst water resources at regional to global scales. This paper presents the first global effort to integrate experimental approaches and large-scale modeling. Using a global soil-moisture monitoring program and a global database of karst spring discharges, the simulations of a preliminary global karstic-groundwater-recharge model are evaluated. It is shown that soil moisture is a crucial variable that better distinguishes recharge dynamics in different climates and for different land cover types. The newly developed dataset of karst spring discharges provides first insights into the wide variability of discharge volumes and recharge areas of different karst springs around the globe. Comparing the model simulations with the newly collected soil-moisture and spring-discharge observations, indicates that (1) improvements of the recharge model are still necessary to obtain a better representation of different land cover types and snow processes, and (2) there is a need to incorporate groundwater dynamics. Applying and strictly evaluating these improvements in the model will finally provide a tool to identify hot spots of current or future water scarcity in the karst regions around the globe, thus supporting national and international water governance.


Introduction
In many countries, karst groundwater is the dominant or even the only available source of freshwater (Stevanović 2019). Climate models indicate that in the next 100 years, karst regions will experience a strong increase of temperature and changes of precipitation in many regions (Hartmann et al. 2014). The potential changes may significantly affect hydrological regimes (Ferguson and Gleeson 2012) and may increase stress on karst water resources. A decrease of water availability can have strong negative impacts on the wellbeing of agriculture, tourism, infrastructure, energy supply, ecosystems and biodiversity. To be prepared, stakeholders and policy makers have to understand the impacts of climate, land use and population change on karst water resources at national and international scales. Policies to ensure an optimal level of adaptation and mitigation can only be developed if quantitative and reliable estimates of potential changes to karst water resources are available at the same scales. Even though strong progress in estimating global water stress was made in the previous years (Wada et al. 2014;Döll et al. 2016;de Graaf et al. 2019), most large-scale modeling studies did not consider the particularities of karst hydrogeology and therefore have limited applicability for water resources management (Hartmann 2016).
The karstic surface and subsurface heterogeneity results in a complex interplay of preferential and diffuse flow patterns. Overall, the hydrological behavior of karst systems shows a duality in its process and storage dynamics (Kiraly 1998): (1) duality of infiltration and recharge processes: diffusive, slow infiltration and recharge into the matrix, and concentrated, Published in the special issue "Five decades of advances in karst hydrogeology".
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10040-020-02258-z) contains supplementary material, which is available to authorized users. rapid infiltration and recharge into the conduits; (2) duality of the subsurface flow field: low flow velocity in the matrix, and fast flow velocity in the karst conduits; (3) duality of discharge conditions-low and continuous discharge during dry periods when the system is dominated by flow through the matrix, and high discharge with high temporal variability during rainfall events when flow through the conduits is dominant. Karstic groundwater flow and discharge have been intensely studied by hydrogeologist (Goldscheider and Drew 2007;Ford and Williams 2013), while recharge generation processes at the shallow subsurface of the karst, i.e. the soil and epikarst, received less attention (Berthelin and Hartmann 2020).
Most karst hydrology models are applied at the scales of individual aquifers (Hartmann et al. 2014) using varying degrees of complexity (Teutsch and Sauter 1991;Sauter et al. 2006;Kovacs and Sauter 2007;Ghasemizadeh et al. 2012;Hartmann et al. 2014). Distributed karst models provide spatially explicit information on groundwater pressure heads and groundwater flow. They are mostly applied at well explored test sites (Chen and Goldscheider 2014;Oehlmann et al. 2014) or were used for theoretical calculations of general behavior of karst hydrology (Covington et al. 2009;Reimann et al. 2014). Lumped karst modeling approaches conceptualize the physical processes at the scale of the whole karst system without being spatially explicit. They consider (1) internal and external runoff (e.g., Jukic and Denic-Jukic 2009;Fiorillo et al. 2015a), (2) epikarst storage and flow processes (e.g., Tritz et al. 2011), (3) groundwater storage and flow in karst conduits and the matrix (e.g. Mazzilli et al. 2019), (4) varying surface and subsurface recharge areas (e.g., Le Moine et al. 2007), and (5) drainage through several springs (e.g., Rimmer and Salingar 2006).
Beyond the scale of individual aquifers, only few studies on quantifying karst water resources can be found. Using observations of specific discharge at multiple sites with high data reliability and precipitation deviations and catchment elevation, Malard et al. (2016) could implement a regional extrapolation of karstic groundwater recharge in Switzerland. Estimating recharge from the difference of mean a n n u a l p r e c i p i t a t i o n a n d m e a n a n n u a l a c t u a l evapotranspiration, Allocca et al. (2014) regionalized karstic groundwater recharge over the southern Apennines in Italy using the areal fractions of limestone and regions without superficial discharge (endorheic areas) as predictors. Huang et al. (2019) showed that terrestrial water storage estimates by the Gravity Recovery and Climate Experiment (GRACE) could be used to quantify the discharge reaction of karst aquifers over the large karst regions of Southwest China.
To predict the impact of climate change and land use changes on karst water availability at larger scales, simulation models are necessary that combine spatial extrapolation or regionalization schemes with the process-oriented model structures. With the aim of quantifying the water balance of the karst dominated island of Crete, Greece, Malagò et al. (2016) developed an extension of the SWAT model (Neitsch et al. 2011) to consider the duality of karstic groundwater. They used a hydrological similarity approach to run their model at the scale of the entire island. Hartmann et al. (2015) used the concept of hydrologic landscapes (Winter 2001) to set up a continental karstic groundwater recharge model over Europe, Northern Africa and the Middle East using a karst specific modelling concept that was previously developed and tested at local scales (Hartmann et al. 2012). Coupled with climate projections (CMIP5, Taylor et al. 2012), the model could be used to estimate future groundwater recharge ).
Yet no approaches to simulate karst water availability exist at the global scale. On the one hand, a lack of observations of karstic groundwater dynamics at the global scale prohibits the extrapolation or regionalization of local information to national or international scales. On the other hand, a lack of conceptual understanding of recharge generation in the karstic shallow subsurface, especially outside the mid-latitude regions of Northern America and Europe, still limits the reliability of large-scale karst recharge models. For those reasons, modeling approaches to provide reliable estimates of karst water resources at the global scale are still not available.
This paper presents the advances of the first global effort to develop a large-scale simulation tool to estimate karst water resources at a global scale to support national and international decision making. Involving large parts of the Commission on Karst Hydrogeology, of the International Association of Hydrogeologists (IAH), an international research project was launched to provide (1) a better understanding of near-surface karst processes by a global soil-moisture monitoring program, (2) new methods to derive regional information karstic of aquifer properties from large numbers of catchment scale observations using a new global database of karst spring discharges, and (3) a systematic approach to incorporate such new understanding into a globally applicable karst simulation model.

Data and methods
Setup of a global monitoring program to characterize soil and epikarst processes Previous work already showed that additional process understanding can be gained by monitoring spatiotemporal variabilities of shallow subsurface hydrodynamics (Penna et al. 2015;Rinderer et al. 2014). Applied in karst regions such approaches can provide more understanding of the local surface heterogeneity and its implication for hydrological modeling. For that reason, a global soil-moisture monitoring program was established to monitor soil-moisture dynamics at a high frequency, at different locations and at different depths. In total, >400 soil-moisture probes were installed across five sites located in Puerto Rico (tropical climate), Spain (Mediterranean climate), the UK (oceanic climate), Germany (mountainous climate), and Australia (semiarid climate). At each site, the probes were split over two different land cover types (forest and grassland) to cover different vegetation cover types.
To account for spatial variability and to minimize the impact of subjectivity when choosing the locations to install the probes, 15 locations for soil profiles were randomly sampled from a uniform distribution at two 20 m × 20 m plots at the forest and the grassland areas of each of the study sites. At each location, vertical profiles with three soil-moisture probes were installed at 5 cm, 10 cm and at the boundary between soil and epikarst (80 cm max). Each profile is connected to a logger that records soil moisture at 15-min resolution ( Fig. 1). In addition, each of the five sites has its own climate station.
Creation of a global database of karst spring discharges to analyses karstic groundwater dynamics Methods to regionalize information from sites with better data availability (see for instance the Prediction in Ungauged Basins (PUB) initiative, Sivapalan 2003;Blöschl et al. 2011) are still limited given the particular complexity of karst systems. Analyzing large data sets of karst system observations would allow for a more comprehensive understanding of regional and global differences of karst system properties. However, an assemblage of karst system observation datasets that would encourage such comparative exercise on larger scales has not been available.
For this reason, the efforts of this study are directed towards the development of a global database of karst spring observations, which would improve access to karst datasets. A framework for the development of the World's Karst Spring (WoKaS) hydrograph database was developed (Fig. 2), involving (1) the identification of karst spring locations, (2) the collection of spring-discharge observations, and (3) the validation of the collected datasets. The previously published "World Karst Aquifer Map" (WOKAM, Chen et al. 2017) was used to support the identification of countries with carbonate rock, karst spring names and locations. An extensive literature review of karst hydrology publications was conducted to further expand the survey range. Discharge observations of the identified karst springs were extracted from publications and national hydrological databases. In addition, a substantial fraction of the observations was provided by individual researchers and members of the IAH Karst Commission. The accuracy and veracity of all collected spring locations, as well as the representativeness of the datasets over the entire globe, were evaluated, which is described in more detail in WoKaS data descriptor (Olarinoye et al. 2020).
For this preliminary analysis, the collected datasets were classified based on elevation, which has been a simple and useful way to compare hydrological system characteristics, especially for analyzing average behavior and variability of recharge and discharge volumes (Stoelzle et al. 2019;Malard et al. 2016). Five classes of springs defined from their elevations in meters above sea level are: L1 ≤ 400 m, 400 m < L2 ≤ 800 m, 800 m < L3 ≤ 1,200 m, 1,200 m < H1 ≤ 1,600 m, and H2 > 1,600 m. The long-term mean discharges and their coefficient of variation (CV) were calculated. Average precipitation values of the spring locations were computed using the GLDAS precipitation datasets (Table 1). With the precipitation information and the simulated recharge values obtained from the model described in the following subsection, the recharge rates and recharge area of those WoKaS springs that had at least 12 months of discharge observations were estimated.

Setup of a preliminary global karst recharge model to quantify water availability
At larger scales the lack of data increases and additional uncertainties arise, since large-scale models are commonly run on grid, while observations are available at point or catchment scale. Hence, a systematic approach to optimize the incorporation of local and catchment-scale karst observations into the development and evaluation of a large-scale karst model is needed. For that reason, a global version of a previously published large-scale karst recharge model (Hartmann et al. 2015) was developed. The model simulates karst recharge processes  Berthelin et al. 2020a, b) based on the general conceptual model of the soil and the epikarst (Fig. 3) accounting for localized runoff, preferential infiltration, evapotranspiration from the soil, and vertical percolation from the epikarst layer towards the groundwater. Actual evapotranspiration (AET) in the model is calculated linearly from the potential evapotranspiration based on the soil saturation. Subtracting AET, the effective precipitation is added to the soil layer. The epikarst receives the vertical flux from soil when it reaches saturation. Groundwater recharge is calculated by a linear relationship dependent on the epikarst storage. When the epikarst reaches saturation, localized runoff towards next model compartment is generated. Groundwater recharge from the epikarst is distributed over diffuse and concentrated recharge using a variable separation factor for conduit and matrix. Details of the model can be found in Hartmann et al. (2013).
In order to incorporate the karstic subsurface heterogeneity, the model assumes the distribution of N compartments to represent the variability of subsurface properties such as soil and epikarst storage capacities, or epikarst hydraulic properties. In the model, these are distributed over N horizontally parallel model compartments (Fig. 3b): S max,i [mm] is the soil or epikarst storage capacity of model compartment i, S max,N [mm] is the overall maximum storage capacity of the soil or the epikarst, K epi,i [days] is the storage constant of the epikarst at model compartment i, K epi,1 [days] is the storage constant of the epikarst at model compartment 1, and a [−] is a dimensionless shape factor. With these equations the water balance of a soil and epikarst layer are calculated at a daily time step in each model compartment. Localized runoff towards model compartments with higher vertical infiltration capacity is initiated when soil and epikarst reach saturation. That way, weak to moderate rainfall events will mostly produce diffuse recharge and/or evapotranspiration, while strong rainfall events will result in concentrated recharge and lower fractions of precipitation are turned into evapotranspiration (Fig. 3b).
Using freely available datasets (Table 1), the model is run over all karst regions in the world (obtaine from Chen et al. 2017;Goldscheider et al. 2020) with daily forcings of precipitation and potential evapotranspiration obtained by the Priestley-Taylor equation (Priestley and Taylor 1972) obtained from (Miralles et al. 2011;Martens et al. 2017). It is run from 1990 to 2019 at a 0.25°× 0.25°spatial resolution where the first 2 years are used as a warm-up period. Other than its application at the continental scale (Hartmann et al. 2015), the preliminary global karst recharge model is not (yet) calibrated with observations of soil moisture and actual evapotranspiration, but it is run with 250 parameter sets sampled from a prior distribution using mean soil and mean epikarst storage capacities of 0-1,250 and 20-700 mm, respectively, mean epikarst storage confidents of 0-50 days and a shape factor a of 0-6. The variability of 250 resulting recharge simulations for each grid cell therefore represents the simulation uncertainty of this preliminary model application.

Evaluation of the global model with the soil moisture and spring-discharge observations
To evaluate the simulated soil storages of the global karst model with the observed soil moisture at the five sites, monthly simulated soil saturation (averaged over the 15 model compartments, Fig. 3b) were compared with the observations at three different depths individually to quantify the strength of their correlation. Since the model simulated soil storage [mm] and the soil-moisture observations provide soil-water content [m 3 /m 3 ], it is not possible to compare them directly without precise knowledge of soil thicknesses and effective porosities. To evaluate the performance of the global recharge model for its soil storage simulations without knowledge about soil thicknesses and effective porosities, instead the simulated and observed soil saturation [−] were compared. The global recharge model can provide the simulated soil saturation by normalizing the simulated soil storage [mm] with the maximum soil storage [mm], which is one of the model parameters. To obtain soil saturation from the observations, the observed soil water content [m 3 /m 3 ] were normalized by the maximum soil-water content [m 3 /m 3 ], which were approximated by the maximum soil moisture in the available records. The mean over all estimated soil saturation time series were calculated for the respective depth class (5 cm, 10 cm, or bottom). Similar to Hartmann et al. (2015) and Sarrazin et al. (2018), the correlation coefficient of the observed and simulated monthly soil saturations was used to evaluate the karst recharge model performance for the three soil depths.
To evaluate the simulated recharge of the global karst model, monthly simulated recharge volumes were compared with mean monthly observed spring discharges of the WoKaS database (Olarinoye et al. 2020, described previously). To minimize the effect of the insufficient length of monthly spring discharge on the correlation, the correlation analysis was only performed for the springs that have at least 12 monthly discharge values (305 springs total). To account for the delay produced by storage and lateral transmission in the phreatic zone, the maximum correlation coefficient of a cross correlation analysis allowing up to 3 months of delay of the observed discharge signal was used to evaluate the simulated recharge. It is assumed that the longer the time delay to the maximum r, the stronger the influence of the phreatic zone.

Results
Over 18 months of soil moisture data were recorded at the sites by the global monitoring program and >400 time series of karst spring discharges were collected for the global karst  (Fig. 4). A 27-year-long record of monthly karstic recharge simulations was produced by the preliminary global model.

Soil-moisture observations at the grassland and forest sites
Starting between April and August 2018, all sites already collected more than 1.5 years of soil-moisture observations. Depending on the site location and the land cover type, they show different patterns in their variability (Fig. 5). The observations at the Puerto Rican site show the highest values of soil moisture correlate to grassland compared to all other sites. However, it is also the site with the lowest soil-moisture values pertaining to forest plots, which is almost similar to the Australian site. The soil moisture is increasing with depth at both vegetation type plots. In particular, at the lowest depth of the forest plot (5 cm), the values of soil moisture are lower in comparison to the Australian site. On the other hand, the deepest probes show values two times higher than the deepest probes in Australia. Considering all depths together, the Australian site shows the lowest soil-moisture values without significant differences between grassland and forest plots. The same is true when considering the soil-moisture variations over different depths.
The soil-moisture variability at the Spanish site is similar to the UK site, however with lower minimum values. The forest and grassland plots do not show significant differences in general, and for all depths considered separately. At the Spanish site, soil moisture tends to increase with depth, which is most visible regarding the grassland plot. At the forest plot, a decrease of average soil moisture is only from 5 to 10 cm, while the soil-moisture variability of the 10 cm and the bottom depth probes are very similar. At the UK site, the soil moisture increases with depth at both sites. The German site shows the highest soil-moisture values after the Puerto Rican grassland plot. At the forest site, soil-moisture values are increasing between the 5 and 10 cm depth and decreasing between 10 cm and the bottom. At the grassland site, the soilmoisture values are decreasing continuously from the surface to the bottom. At both the German grassland and the forest sites, the deepest probes show the largest spread in their soilmoisture dynamics. More information about the temporal dynamics of soil moisture at the different sites is provided in Fig.  S1 of the electronic supplementary material (ESM).

Collected karst spring hydrographs data
Through the established data collection framework and a combined community-effort, the WoKaS database presently archives more than 400 karst spring-discharge observations globally (Olarinoye et al. 2020). The length of the datasets ranges from a few months up to 120 years with a median record length of 14 years (Table 2). In all, 50% of the datasets contain discharge records sampled at a daily or subdaily frequency but datasets in upper quartile have an observation temporal resolution of 4 days and above, most of which are datasets with longer data records. On average, 95% of the datasets in the WoKaS database provide continuous discharge records.
The average discharge of collected karst springs for the five elevation classes spreads across 10 −4 to 10 2 orders of magnitude (Fig. 6). Larger springs are located at lower elevations up to 1,200 m. (elevation class L1, L2 and L3). Most springs ) have lower discharges. Springs located at lower elevations (L1, L2 and L3) show higher CVs compared to those located at higher elevations (H1 and H2) with less variability among different springs. From the recharge model described in subsection 'Setup of a preliminary global karst recharge model to quantifywater availability' recharge values from approximately 300 spring locations were obtained. Therefore, the recharge rate and recharge area analyses (Fig. 6c, d) are provided for the subset of WoKaS datasets for which recharge values have been estimated. The recharge rates range from low to very high values. A systematic pattern between recharge rates and elevation is found. High recharge rates of up to 70% are observed among L1, L2 and L3 springs (Fig. 6c). In all, 50% of the low-elevation springs (L1-L3) have a recharge rate higher than 45%, while the high-elevation springs (H1-H2) within the same quantile have recharge rates >30%. Irrespective of the elevation, the estimated values show a high variability in the recharge rates. In Fig. 6d, extreme ranges of recharge areas from values <1 km 2 to larger areas of up to 10 4 km 2 are shown. Unlike for the recharge rates, there is no systematic pattern or order found between the recharge area and elevation; however, all spring classes have an almost similar range of median values which is slightly less than 100-km 2 recharge area. About a quarter or slightly more of the recharge areas at all classes are <10 km 2 , and at least the upper quartile or even more have areas >100 km 2 .

Global groundwater recharge simulations
The mean annual recharge volumes derived for the period 1992-2019 resemble the meteoric water availability in the different regions in the world (Fig. 7a). Rainy regions such as Scotland and Ireland, coastal regions and monsoonal regions are also characterized by recharge volumes close to 1,000 mm/year or more. On the other hand, regions that are characterized by aridity show average recharge volumes as low as just few mm per year such as in Northern Africa, central Northern America, the Middle East or the Himalayas. In the same regions, model uncertainty tends to larger values, with standard deviations as large or even larger than the average annual recharge (Fig. 7b), while uncertainty remains low in the wetter regions.

Evaluation of the global model with the soil moisture and spring-discharge observations
The simulations of soil saturation of the global karst recharge model are compared with the observed soil-moisture dynamics at the five sites. At its present state, the model tends to overestimate the monthly average soil saturation at Austrian,  German and UK sites regardless of the land types (Fig. 8). For Puerto Rico, the soil saturation of forest is overestimated, as well. Generally, linear relationships with varying slopes between observed and simulated monthly average soil saturation   (Table 3). In addition, soil saturation shows different variability. Especially at Puerto Rican sites, it spreads in different soil saturation ranges with forest <0.4 and grassland between 0.4 and 0.8 (Fig. 8).
Overall, the correlation coefficients r of the monthly observed and simulated soil saturation reach values up to 0.94. Weak relationships, r < 0.45, go along with insignificant correlation (Table 3). Forest and grassland show different strengths of correlation, with stronger correlation for the forest as compared with that of grassland (except for the UK site).
The correlations between the monthly observed karst spring discharge and the corresponding monthly simulated recharge ( Fig. 9) show that 47 and 59% of the springs show a correlation coefficients r ≥ 0.5 without and with consideration of time delay from recharge to discharge, respectively. The larger the value of r, the more significant the correlation is observed. As expected, the correlation between recharge and discharge can represent how strong the recharge is linked to the discharge. It can also be seen that the time delay from recharge to discharge helps to obtain a better correlation for some regions (Fig. 9). Few negative correlations between recharge and discharge suggest that local conditions of springs, e.g. the topography, could substantially affect this relationship.

A better characterization of karstic recharge processes by soil-moisture dynamics
The dynamics of the collected soil-moisture observations allow for preliminary interpretation and new region-specific and  land use-specific insights (Fig. 5). A strong linkage of climate and soil moisture is found. For instance, the highest soilmoisture values occur at the site with a tropical climate at the grassland plot reflecting the wet tropical climate conditions. At the forest plot, rather low soil-moisture values can be explained by the dense network of tree roots and small amount of soil that can store the infiltrating water. High values of soil moisture are also measured at the German site, where high annual volumes of precipitation prevail. On the other hand, the low soil-moisture values measured at the Australian site are coherent with the semiarid local climate conditions. Despite their different climatic regions, the Spanish site (Mediterranean climate) and the UK site (oceanic climate) show similar variability of observed soil-moisture dynamics. This is probably due to their similar annual precipitation volumes of 760 and 815 mm/year for the Spanish and UK site, respectively, and to their mean annual temperatures of 14°C (ES) and 5.4-14°C (UK) yet occurring with different strength of seasonality (Berthelin et al. 2020b). The control of climate on soil-moisture dynamics, and vice-versa, is well-known (Seneviratne et al. 2010) but in order to derive improved concepts of groundwater recharge processes from soil-moisture dynamics, more parameters have to be considered such as soil texture, antecedent moisture conditions, vegetation, and the epikarst (e.g., Perrin et al. 2003;Heilman et al. 2014;Fu et al. 2015;Martos-Rosillo et al. 2015). Comparing the evolution of soil moisture with depth, the probes at 5-cm depth present the lowest values at every site, and soil moisture is increasing with depth. This is most probably linked to evaporation processes that have a stronger impact on shallow soil-water storage (Martini et al. 2015;Sprenger et al. 2016). Only the German site presents soil-moisture values that decrease with depth indicating rapid shallow subsurface flow paths (Chifflard et al. 2019), which may be favored by the strong slopes of this site and its location in the mountains.
Yet, the comparison between sites, different soil depths and land cover types remains qualitative and preliminary. The three main parameters explored in the preceding (climate, land cover and depth) are not the only ones that influence soilmoisture dynamics. In addition, the climate could affect soilmoisture dynamics differently in different seasons (Berthelin et al. 2020b) and might be dependent on precipitation amount as well as intensities. The influence of antecedent soilmoisture conditions on recharge initiation could be revealed by considering a larger number of extracted soil-moisture events and their pre-event soil storages (Demand et al. 2019). At those sites, where observations of groundwater, or of related fluxes like stream, discharge, spring discharge or drip in caves are available, methods to estimate recharge from soil-moisture observations by simple models (Baker et al. 2020) or data-driven approaches can be explored (Arnold et al. 2020). Those approaches may be supported by analysis of stable isotopes in soil water as already proven to be useful in non-karstic settings by Sprenger et al. (2015). Overall, with another 18-24 months of monitoring at the five sites, a dataset will be provided to advance the conceptual understanding of karstic recharge and evapotranspiration processed both qualitatively and quantitatively.

Pathways to upscale local understanding by the WoKaS database
The WoKaS database tends to contain larger springs located at lower elevations (Fig. 6). Hydrologically, springs at lower elevations are located at or close to catchment outlet; therefore, they drain a larger catchment area producing the large discharge volumes (Kresic and Stevanovic 2009). Similarly, a higher and wider range of CV values is associated with spring discharges at lower elevations. This implies that springs at higher elevation have more consistent discharge variability throughout the data record period, which may be due to the seasonality produced by snow accumulation and snow melt (Chen et al. 2018). Since springs at lower elevation drain a larger catchment area, the recharge area is consequently large with variable recharge sources. This and other climate variables could be attributed to the higher discharge variability of springs at lower elevation.
The high recharge rates up to 70% found at WoKaS springs' locations is no surprise. Groundwater recharge is known to be higher in karst areas compared to other landscapes  where more large fractions of the total precipitation volume can infiltrate into groundwater (Bonacci 2001;Fiorillo et al. 2015b). Usually, higher elevations receive more precipitation and higher recharge rates would be expected as Fig. 9 Distributions of the correlation coefficients r between the monthly simulated recharge and monthly observed karst spring discharge (from WoKas, Olarinoye et al. 2020). Blue and orange bars represent the correlation without and with time delay from recharge to discharge, respectively well. This was found, e.g., in the Swiss Alps by Malard et al. (2016) or the Italian Apennines by Allocca et al. (2014); however, an increase of recharge rates with elevations does not occur in the global dataset as it also covers mountain ranges in very dry climate regions such as Central Northern America, the Middle East and Southern Australia (Fig. 7a). Considering the range of the corresponding recharge areas (obtained by water balance, see section 'Creation of a global database of karst spring discharges to analyses karstic groundwater dynamics'), similar variability and averages for all elevations is found, indicating that the dataset is not biased towards different scales of karst systems at different elevations.
The present analysis only gives an overview of the attributes and characteristics of karst springs by exploring the collected datasets. The database still provides lots of potentials yet to be explored. In future analysis, the dynamics of karst springs in different regions will be explored to see how local factors influence discharge and recharge variability. The expected outcome of this analysis will enable us to identify important local drivers and even predict spring behavior in regions with nonreliable or no observation records. Also, the estimated recharge areas could be a first step for their spatially explicit delineation (Malard et al. 2015). As springs also reflect the dynamic behavior of karst aquifers, important information such as recession parameters derived from the large datasets could be used to infer the dominance of conduit and matrix contributions in different regions. Presently, the WoKaS datasets are is available for download in an online stationary repository, but efforts will be made to provide the datasets directly through a web platform with a graphical user interface. Such development will allow for continuous growing of the database, adding other complementing datasets and a web tool for instant analysis.

Model deficiencies revealed by evaluation with the newly collected observations
The simulated mean annual recharge volumes mostly reflect the regional climatic conditions (Fig. 7a), a result which is very similar to its previous continental-scale application over Europe Northern Africa and the Middle East (Hartmann et al. 2015). Small differences in simulated average recharge volumes are most probably due to a new delineation of karst areas-global model: WOKAM, Chen et al. (2017); continental model: "Global distribution of carbonate rocks", Williams and Ford (2006)-and different simulated time periods (global model: 1992-2019continental model: 2002. However, when looking at the simulation uncertainties (Fig.  7b), the preliminary character of the global model is more obvious. Especially in arid regions, the simulation uncertainty exceeds 100% making simulations of karstic groundwater recharge basically useless for water management in those regions. Yet, simulation uncertainty strongly reduces in semiarid wetter regions where even these preliminary simulations could be useful for water managers and water governance. In those regions, the fractions of precipitation turned into recharge are substantially higher compared to arid regions making precipitation itself a good predictor of groundwater recharge and reducing the relative impact of the uncertain preliminary model on the mean annual recharge estimates.
Through the comparison between observed and simulated soil saturation (Fig. 8), an obvious deviation of simulations of the global model and the observations can be seen, mostly expressed through an overestimation of soil saturation by the model. This deviation is influenced by several aspects. The simulated soil saturation is averaged over a large grid that represents the integral response for this large area, while the observed soil saturation is measured at a specific point that can differ a lot because of heterogeneities of soil properties and land cover from site to site, i.e., there is a problem of incommensurability (Beven 2018). Considering the coefficient of correlation between simulations and observations as a measure of model performance (similar to Hartmann et al. 2015;Sarrazin et al. 2018), this problem is partially circumvented as r is not affected by differences of effective porosities. Comparing the coefficients of correlation for the different sites and different land cover types (Table 3), it can clearly be seen that the model performs well for the Spanish forest and grassland sites and the Puerto Rican forest site. Bad correlations that are sometimes both even significant, are found at the Australian and German sites for both land covers, and the Puerto Rican grassland site. The different performances between grassland and forest point towards the very simplified representation of land cover in this preliminary model (Sarrazin et al. 2018). While the weak performance at the German site, which is located at~1,450 m above sea level, is most probably due to the neglecting of snow processes in the model, the model deficiencies at the Austrian site could be due to general uncertainty of the gridded input data for this region as already discussed by Baker et al. (2020).
Considering the correlation between simulated monthly recharge and observed WoKaS spring discharge, a large number of relatively high r values are found, despite the preliminary state of the model (Fig. 9). However, there are also a substantial number of springs with weak linear relationships and even negative correlations between simulated recharge and observed discharge. This could be explained by the limited consideration of the location and size of the recharge area in the model. Since it is not possible to delineate the real recharge area of every spring, the simulated recharge of the grid cell where the spring is located is used as the recharge of this spring. However, the recharge area ranges across several grid cells, which may differ strongly from its topographic area (Le Moine et al. 2007;Longenecker et al. 2017;Le Mesnil et al. 2020). Due to this difference, the correlation for these springs can be biased. Another, even more probable reason for the weak correlations is the lack of groundwater processes in the preliminary model. This is confirmed by the improved correlations between recharge and discharge that is obtained after allowing for the time delay from recharge to discharge.

Towards reliable simulations for (inter) national water governance in karst regions
The comparison of simulated and observed soil moisture clearly indicates that land cover has significant influence on soil moisture as well as evapotranspiration mentioned already. Land cover affects the partitioning of precipitation into evapotranspiration, soil moisture, and surface runoff. This highlights the importance of including explicit land use types to improve global karst recharge modelling, allowing to investigate impacts of land use change on the recharge and discharge (Sarrazin et al. 2018). The poor performance at the mountain site in Germany shows the need to add a snow model in order to include karst regions located in mountain regions (Chen et al. 2018). More recent global input datasets such as MSWEP V2 (Beck et al. 2019) will help to improve the recharge simulations at dry sites such as the Australian site. A need to include a karstic groundwater model is revealed through introducing a time delay between recharge and discharge (Fig. 9). The improved correlation between simulated recharge and observed discharge after introducing such delay suggests that, despite the fast karstic flow paths, also slow groundwater transmission and storage takes place in the phreatic zone. Adding a groundwater routine that considers system properties such as the distribution of the conduit networks, and the permeability of the matrix, will provide a better representation of the delayed response of karst springs to a recharge signal (Geyer et al. 2008;Covington et al. 2009).
The recharge area of karst aquifers is the most common spatial unit to investigate and model karst springs; however, larger river basins that drain karst regions are often partially covered by nonkarstic areas. Water management at these basins therefore needs to understand the combined behavior of both systems. Only few studies (e.g., Rimmer and Salingar 2006;Chen et al. 2018) have considered both karstic and nonkarstic components in catchment-scale modeling. Challenges remain for modeling such systems such as intercatchment groundwater flow can cross the topographic boundary of a catchment and result in unclosed water balances (Le Mesnil et al. 2020). Neglecting this disagreement of surface and subsurface catchments will limit the representation of karstic and nonkarstic hydrologic processes in combined modelling systems. Therefore, identification and quantification of inter-catchment groundwater flow is of great importance. This may be achieved by diagnostic signatures based on independent datasets and water balance (e.g., Liu et al. 2020) and new approaches to integrate this information into regional models with combined karstic and nonkarstic processes representations.

Conclusions
This paper shows the most recent advances in developing a global karst modeling system using a global soil-moisture monitoring program and a global database of karst spring hydrographs. Comparing the simulations of a preliminary version of the first global karst recharge model with the soil-moisture observations reveals that improvements of the soil and epikarst processed in the model are still necessary to obtain a better representation of different land cover types and snow processes. The comparison of observed spring discharge with the simulated recharge values strongly points towards the need to incorporate groundwater dynamics including the interplay of partially overlapping surface and subsurface catchments and the influence of nonkarstic units in karst-dominated river basins. Consequently, the comparison of the preliminary model with the newly collected soil-moisture data and spring-discharge observations provides detailed and explicit directions to make important advancements towards the first global karst simulation model. Such modeling system will not only provide information about water availability in the simulated catchments. Karst aquifers provide drinking water for a large part of the world population (Ford and Williams 2013) and are among those groundwater resources that are far from being overexploited (Stevanović 2019). Applied at a global scale and fed by climate projections, the model will also allow to identify hot spots of current or future water scarcity in the karst regions around the globe and where karst aquifers may mitigate water shortages. That way, it can support national to international water governance to develop regional and local mitigation measures to successfully tackle the impacts of climate change, land use a change and population growth. otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.