Global‐Scale Groundwater Recharge Modeling Is Improved by Tuning Against Ground‐Based Estimates for Karst and Non‐Karst Areas

Quantification of groundwater recharge is highly important for understanding freshwater systems and sustainable management of both groundwater and surface water bodies, but it is very uncertain. To improve the global‐scale simulation of diffuse groundwater recharge (GWR), we developed, for the first time in global hydrological modeling, an approach for explicitly simulating GWR in karst areas, based on the World Karst Aquifer Map and ground‐based estimates of karst GWR. Ground‐based estimates of both karst and non‐karst GWR from 5,600 locations were aggregated to produce estimates for 816 globally distributed 0.5° grid cells, including 64 representing karst GWR. Applying WaterGAP2.2e, we found that karst GWR is best approximated by total runoff from land (without soil overflow and urban runoff). The GWR algorithm for non‐karst GWR was tuned against ground‐based estimates for 422 non‐karst grid cells outside of Australia. While simulated GWR tends to underestimate GWR in non‐karst areas with mean annual ground‐based GWR > 200 mm/yr, it overestimates national GWR in six out of seven European countries. With an increase in the coefficient for the discharge from groundwater to surface water bodies, the updated GWR algorithm results in an improved fit of simulated streamflow to observations, including low flows. Global mean annual GWR and thus renewable groundwater resources (without Greenland and Antarctica) are estimated to be about one‐half of the total renewable water resources, or 21,000 km3/yr, which is over 50% higher than the standard model estimate. 20% of global GWR is generated in karst areas, which cover about 11% of the world's land surface.


Introduction
Groundwater is a crucial source of water supply for humanity, accounting for one-quarter of total freshwater abstractions worldwide in 2017, and even more than one-third of the abstractions for domestic use (IGRAC, 2018;Müller Schmied et al., 2023).Groundwater abstractions lead to a decline in the groundwater table and a decrease in groundwater discharge to streams (de Graaf et al., 2019), the latter causing reduced seasonal low flows (Wittenberg, 2015).The decline of the groundwater table and discharge negatively impacts freshwater biota (Saccò et al., 2024) and increases the ratio of groundwater abstractions to groundwater recharge (GWR).Renewable groundwater use is a dynamically stable re-equilibrium of abstraction and GWR in the long term (Cuthbert et al., 2023).Accurate quantification of GWR is thus an essential foundation for sustainable groundwater management.
GWR is defined as the addition of water to the water-saturated zone below the groundwater table, either from the overlying unsaturated zone, known as diffuse GWR, or from a surface water body, termed as point GWR (Healy, 2010).Influenced by the complex nature of soil, terrain, subsurface rock structure, groundwater-surface water hyporheic exchange, and climate variation, GWR is known to be highly heterogeneous in both space and time (Fu et al., 2019).While streamflow measurements provide reliable estimates of the total renewable water resources in the upstream area, GWR cannot be directly measured.Indirect methods of GWR measurement, including sediment profile analysis, water table fluctuation, hydrograph separation, etc., yield mostly local and uncertain estimates (Levin et al., 2023).These observation uncertainties arise not only from measurement errors but also from the inherent complexity of the hydrogeological environment.We hereafter refer to these indirect GWR measurements as "ground-based GWR estimates."Continental-to global-scale GWR assessments are mainly performed using global hydrological models (GHM) (Döll & Fiedler, 2008;Reinecke et al., 2021).Most GHMs only simulate diffuse GWR from soil zones (e.g., Burek et al., 2019, for CWatM;Verkaik et al., 2022, for PCR-GLOBWB), while GWR from surface water bodies such as rivers, lakes, and wetlands, which depends on the hydraulic gradient between surface water and groundwater, is currently not simulated by most GHMs.In the following, we only consider diffuse GWR and refer to it simply as GWR.Döll and Fiedler (2008) were the first to estimate global-scale GWR by applying the GHM WaterGAP, where daily GWR is computed as a fraction of daily total runoff.The fraction is calculated as a function of relief (based on slope), soil texture, hydrogeology, and the occurrence of permafrost and glaciers as well as aridity.Out of the 16 GHMs whose model structure was compared in Telteu et al. (2021), two others (MATSIRO and WAYS) follow the GWR approach of WaterGAP.Most of the other GHMs first simulate surface runoff and/or interflow before computing GWR as a function of the saturation of the lowest soil layer or use matric potential and soil hydraulic conductivity to simulate the location of and the flow across the groundwater table, before simulating surface runoff and interflow (Telteu et al., 2021, their Table S30).Some take into account capillary rise.The simulated GWR varies significantly among the GHMs (Reinecke et al., 2021;West et al., 2023) and depends also on the applied climate forcing (Döll & Fiedler, 2008).GWR estimates of GHMs have rarely been compared to ground-based estimates.Less than 60 globally distributed estimates were used for evaluating GWR computed by WaterGAP (Döll & Fiedler, 2008) and PCR-GLOBWB (Jasechko et al., 2014).West et al. (2023) compared 124 ground-based estimates in Africa compiled by MacDonald et al. (2021) to the output of 8 GHMs.They found that among the 124 × 8 comparisons between simulated GWR and ground-based estimates, GHMs underestimated ground-based estimates in 56% of the cases and overestimated them in 18% of the cases, while 26% fell within the uncertainty range of the ground-based estimates.GWR from H08, MATSIRO, and WaterGAP showed the lowest discrepancies, and WaterGAP tended to underestimate GWR in all four investigated climatic regions.
GHMs serve to estimate GWR dynamics at the sub-annual to decadal timescales, and also under the impact of climate change.If the focus is only on the recent long-term GWR and their spatial distribution, it is possible to upscale local ground-based GWR estimates to the global or continental land area.This upscaling is achieved through statistical modeling that incorporates spatial data on relevant predictors, as demonstrated in global studies (Berghuijs et al., 2022;Mohan et al., 2018) and research focused on Africa (MacDonald et al., 2021).However, this upscaling approach suffers from inconsistent ground-based GWR estimates in temporal scope and geographic coverage, as well as the lack of data in extensive areas across the globe (Moeck et al., 2020).For example, Mohan et al. (2018) were suspected to overestimate GWR at high northern latitudes as no GWR data for permafrost areas were available (Reinecke et al., 2021).The upscaling studies found that spatially distributed GWR is best predicted by long-term averages of either precipitation only (MacDonald et al., 2021), aridity only (Berghuijs et al., 2022) or precipitation, potential evapotranspiration and land cover (Mohan et al., 2018), while other characteristics that are known to affect GWR locally, for example, relief and soil characteristics, were not identified as relevant predictors.Additionally, none of the GHMs or upscaling studies considers the particularities of karst hydrology, which includes substantially high infiltration and GWR compared to non-karst terrains (Ford & Williams, 2007;Stevanović, 2015).
The main reason for neglecting karst in GHMs was that global geological maps did not indicate karst.While still not indicating karst but karstifiable areas, the World Karst Aquifer Map (WOKAM) (BGR et al., 2017;Chen et al., 2017) can now serve as a basis for taking into account karst in GWR modeling.However, GWR observations in karst are rare (Hartmann et al., 2015) and were not included in the Moeck et al. (2020) data set.While the World Karst Spring hydrograph database (Olarinoye et al., 2020) provides a compilation of >400 time series of karst spring discharge mainly in Europe and the US, the lack of information on the spring-specific recharge areas prevents the derivation of GWR in volume per area, which is necessary for comparison with other areas (Ford & Williams, 2007).Carbonate and, less importantly globally, evaporite rocks are karstifiable, for example, large conduits can form by the dissolution of solid rock into the flowing water, both at the surface and in the subsurface.The degree of karstification of carbonate rock varies strongly.GWR is generated by the movement of water through the epikarst, both as diffuse recharge like in non-karstified areas and as concentrated recharge from lateral movement of surface and subsurface water into karst conduits, for example, into dolines formed by carbonate dissolution (Hartmann et al., 2014).The conduit flow passes the root zone quickly, thus avoiding evapotranspiration.Karst aquifers contribute to water supply in more than 140 countries, even though their flashy discharge regime and the related high vulnerability to pollution are problematic, and it was estimated that they are the source of ca.13% of global groundwater abstractions (Stevanović, 2019).Neglecting the specifics of karst GWR results in systematically underestimating GWR in karst areas, for example, in the Mediterranean (Hartmann et al., 2013(Hartmann et al., , 2015)).
Hydrological modeling of karst requires local observations for parameterization, and Hartmann et al. (2013) were not successful in transferring model parameters to ungauged karst aquifers.To quantify karst GWR at the continental scale, Hartmann et al. (2012Hartmann et al. ( , 2015) ) developed the VarKarst-R model, where the sub-grid spatial variability of the soil-epikarst compartment is taken into account and lateral flow within the grid cell can lead to concentrated conduit flow.For their study in Europe and the Mediterranean, with a spatial resolution of 0.25 by 0.25°, they initially used a vast ensemble of 25,000 parameter sets, refined by comparing with observations of actual evapotranspiration from FLUXNET (Baldocchi et al., 2001) and soil moisture from ISMN (Dorigo et al., 2011), to select the best 250 sets for assessing karst GWR.However, such complexity of selecting locationspecific parameter sets, precise parameterization, and the requirement for extensive local observations make this approach less practical for GHMs, which require more generalized and scalable methodologies.Additionally, the coverage of karst landscapes outside of Europe with FLUXNET stations is limited.
In this paper, we present how karst is taken into account when modeling GWR in the GHM WaterGAP.In addition, we describe how the computation of GWR in non-karst regions was adjusted using up-to-date databases of land surface slope and soil properties.Both algorithm modifications were derived by optimizing the fit of simulated GWR values to a large number of ground-based estimates assigned to 0.5°grid cells.Finally, we analyze how the methodological advances of GWR modeling impact the estimate of global long-term average GWR, its spatial distribution, and the contribution of karst to total GWR, and how they affect the fit of simulated streamflow to observations.In the last two sections, we discuss the GWR computation and model outputs and draw conclusions.

Methods and Data
This study relies on a compilation of globally distributed ground-based GWR estimates, sourced from existing compilations that include and distinguish estimates of karst GWR, with an assignment of best ground-based estimates to 0.5°grid cells (Section 2.1).These independent ground-based GWR estimates are utilized in conjunction with the output of the GHM WaterGAP 2.2e (Section 2.2), particularly runoff from land, and new geographic data sets including WOKAM to optimize the computation of GWR (Sections 2.3-2.7).

Collating Ground-Based GWR Estimates
Ground-based GWR estimates can be derived from measurements in the saturated and unsaturated zones and in streams.Estimates based on measurements in the saturated zone and in streams include the effect of capillary rise, while measurements in the unsaturated zone quantify potential GWR that may or may not reach the groundwater table .To obtain reliable GWR estimates, Scanlon et al. (2002) recommended combining several GWR measurement methods.Local ground-based GWR estimation methods encompass chloride mass balance, tracer analysis (e.g., chlorine and tritium), the analysis of groundwater table fluctuation, and various other methods (e.g., MacDonald et al., 2021;Moeck et al., 2020).The first two types of methods were used to determine about 80% of the >5,000 local GWR estimates compiled by Moeck et al. (2020).GWR at the scale of drainage basins can be estimated by stream hydrograph separation and by the water budget method (Huang et al., 2023;Levin et al., 2023).Some ground-based estimation methods include hydrological modeling (e.g., Eilers et al., 2007).
We compiled a global-scale data set of ground-based GWR estimates by synthesizing six major data sets and by identifying additional ground-based estimates for karst from individual publications (summarized in Table 1).To enable a detailed model evaluation in karst regions, we reviewed each data source in the data set of Hartmann et al. (2015) by reviewing each data source.Three out of the originally 36 estimates (for Afka in Lebanon, Aydincik in Turkey, and Sanità in Italy) were eliminated due to their very high value compared to WaterGAP precipitation (GWR larger than model precipitation).To obtain more ground-based karst GWR estimates, a literature search and review was conducted using Google Scholar and CNKI net (for Chinese studies) with the keywords "karst," "carbonate," "dolomite," "limestone," and "chalk."This resulted in an additional 31 groundbased karst GWR estimates.
When generating the new compilation of ground-based GWR estimates, duplicates were first removed by comparing the location, reference source, and data set comments.To ensure all estimates represent natural diffuse GWR, we excluded study sites affected by irrigation, in urban areas and beneath surface water bodies (ephemeral or perennial streams and lakes); this step was only applied to the Scanlon et al. (2006) data set as other data sets have not provided such information.To prevent GWR discrepancies caused by differences between the precipitation indicated for the ground-based GWR estimates and the precipitation driving the WaterGAP model, only estimates with a recharge coefficient (ratio of GWR to mean annual WaterGAP precipitation) of less than 80% were considered comparable, and thus retained.This removal of observation is only applied to non-karst areas.
The ground-based GWR estimates were allocated to the 0.5°by 0.5°grid cells of WaterGAP based on the latitude and longitude information in the data source.In the case where the deviation between locally recorded precipitation and WaterGAP precipitation input exceeds ±50%, the geolocations of these estimates were checked using Google Maps, and if necessary corrected.The "best" ground-based estimates for half-degree grid cells were determined by applying the following priorities: (a) estimates of recharge to karst aquifer (only 2 cells had both karst and non-karst observations, and the non-karst observations were neglected in this study); (b) estimates at a spatial scale that is closest to the grid size (i.e., ∼2,500 km 2 ); (c) long-term (multi-decadal) estimates take precedence over shorter period estimates; (d) estimates using chloride mass balance (CMB) method than other methods such as environmental tracers, water-table fluctuation, as CMB offers useful insights of recharge over the long-term (only applied for Scanlon et al., 2006 data set).If a best estimate could not be identified, we adopted the arithmetic mean of all records distributed in the same grid cell as the ground-based estimate.The Supporting Information S1 (e.g., mean precipitation and temperature, soil group, bulk density, aridity, recharge estimation method(s), study scale, and reference source) was retained whenever possible.
In summary, we identified, from 341 individual studies, a total of 816 ground-based estimates of GWR globally that are assumed to be representative for 0.5°grid cells (Figure 1).Sixty-four of these estimates are from karst studies but with 3 estimates located in cells that are not identified as karst aquifers according to the World Karst Aquifer Mapping (see Section 3.2).The minimum and maximum GWR values, either reflecting spatial variability, assumptions made, or measurement methods applied, were identified for 545 out of the 816 0.5°cells, but only the grid-level best estimates were used for the final comparison.

Simulating GWR With WaterGAP 2.2e
WaterGAP version 2.2 is a GHM that simulates freshwater fluxes and storages as well as human water use from both groundwater and surface water bodies on all continents of the Earth except Antarctica, in 67,420 0.5°grid cells, and with a daily time step.In this study, WaterGAP 2.2e was run with the global daily GSWP-W5E5 climate forcing (Lange et al., 2022), which is a merge of the GSWP3 data set (Dirmeyer et al., 2006) for the period 1901-1978 and the W5E5 data set (Lange et al., 2021) for 1979-2019.WaterGAP simulates the hydrological processes on continents by a sequence of water balances for 10 storage compartments, as depicted in Figure 2, where the processes influencing (diffuse) GWR simulation are illustrated.The vertical water balance of each grid cell includes the canopy, snow, and soil water balances.The lateral water balance includes the groundwater balance within each grid cell and the balance of several surface water types among the cells (lakes, wetlands, man-made reservoirs, and rivers).Human water use is assumed to affect only the water storages in the lateral water balance.WaterGAP 2.2e is calibrated against the mean annual streamflow by adjusting model parameters in the upstream grid cells of 1,509 globally distributed streamflow gauging stations.This assures that, for most calibration basins, the long-term total runoff of each calibrated basin is within 10% of the observed value (Müller Schmied et al., 2021).
Before the computation of (diffuse) GWR, the total runoff from land R total is computed, that is, the part of the precipitation and snowmelt that reaches the soil (P eff ) after the computation of the canopy and snow water balances that is neither stored nor evapotranspirated.Total runoff from land has three components.On the urban fraction of each grid cell, 50% of P eff is directly turned into urban runoff R 1 .Then, soil water overflow R 2 is computed; if the sum of soil water storage S s and P eff exceeds S s,max , the cell-specific maximum soil water content (mm) which depends on landcover and soil type, the exceeding amount becomes R 2 .Finally, the runoff component R 3 (mm/d) is calculated following Bergström (1995) with where γ is the runoff coefficient ( ), which is adjusted by model calibration.
Daily GWR R g is estimated as part of R 3 following a heuristic approach based on expert knowledge of the relation between GWR and the spatially varying physical characteristics relief, soil texture, hydrogeology, and the existence of glaciers or permafrost.Capillary rise is not simulated.Daily GWR is computed as a fraction of R 3 and is constrained by a maximum daily GWR related to soil texture, with where R gmax is the soil texture-specific maximum groundwater recharge rate, with values of 7, 4.5, and 2.5 mm/d for coarse-, medium-, and fine-textured soils, respectively.f g is the groundwater recharge factor ( ) ranging between 0 and 1, with where f r , f t , f h , and f pg are relief-, soil texture-, hydrogeology-, and permafrost/glacier-related factors ( ), each in the range between 0 and 1.The values of the four factors were related to the grid cell-specific characteristics in the form of look-up tables, using expert guesses.With steeper slopes, finer soil textures, and less permeable aquifers, f g is expected to decrease, while glaciers and permafrost are assumed to prevent groundwater recharge (Döll & Fiedler, 2008).The heuristic approach in WaterGAP reflects the impact factors that can be seen in the large GWR data set from Moeck et al. (2020).They found decreased GWR with increasing relief.The sum of R 1 , R 2 and the fraction of R 3 that does not recharge the groundwater becomes fast surface runoff.The simulated R g is added to the groundwater compartment and can leave the grid cell as discharge to the surface water bodies in the grid cell or as net groundwater abstraction.Groundwater discharge is computed as the product of groundwater storage and a globally constant groundwater discharge coefficient (Figure 2).WaterGAP further computes focused recharge beneath surface water bodies in (semi)arid grid cells, which is not considered in this study.
Based on additional comparison between groundwater simulation and independent estimates of GWR and groundwater depletion performed in Döll and Fiedler (2008) and Döll et al. (2014), respectively, regional adaptations have been made.(a) If a grid cell is (semi)arid (see Figure 1 for global distribution) and has a medium to coarse texture, it is assumed that daily recharge R g occurs only when daily precipitation exceeds a critical value of 12.5 mm/d ("(semi)arid-coarse adjustment"), also to implement the observation groundwater recharge in dry areas only happens in case of heavy rainfall; (b) For the Mississippi Embayment, the factor f g was reduced from 0.8 to 0.9 to 0.1 as groundwater depletion was strongly underestimated as compared to in-situ estimates, suggesting overestimated GWR (Döll et al., 2014); (c) The wetland that encompassed all of Bangladesh in the original wetland data set used in WaterGAP was removed as in WaterGAP there is no groundwater recharge under wetlands in humid areas.
WaterGAP uses three different types of areas, that is, cell area (time-invariant), continental area (time-invariant), and land area (time-variant).The cell area is the area of 0.5°grid cells calculated with an equal area cylindrical projection.The continental area is the cell area minus the ocean area.The land area is the continental area minus the area of all surface water bodies except the rivers.The land area in each cell varies daily as the area of the surface water bodies is adjusted with the daily varying water storage in these water bodies.All computations of the vertical water balance are done for the land area of each grid cell.So GWR is computed, in units of mm/d, for the land areas, but converted to the continental area as standard WaterGAP output (resulting in lower values in grid cells with surface water bodies).To ensure comparability to ground-based GWR estimates, all WaterGAP results in this paper are presented with respect to land area.Water Resources Research 10.1029/2023WR036182 WAN ET AL.

Localization of Karst Areas
We determined the areal fraction of karst f karst ( ) in each 0.5°grid cell as the fraction of continental area with exposed carbonate and evaporite rocks, by using the information from the World Karst Aquifer Map WOKAM (BGR et al., 2017;Chen et al., 2017).While the map only provides polygons that indicate karstifiable areas but not karst, it is reasonable to assume that most exposed carbonate and evaporite rocks are karstified (Chen et al., 2017).In WOKAM, five categories of karstifiable rock are distinguished: discontinuous carbonate rocks, discontinuous evaporite rocks, continuous carbonate rocks, continuous evaporite rocks, and mixed carbonate and evaporite rocks.Areas with more than 65% of carbonate or evaporite rocks were mapped as "continuous" and areas with 15%-65% as "discontinuous."Areas that contain more than 15% of each rock type were mapped as "mixed carbonate and evaporite rocks" (BGR et al., 2017).For each 0.5°grid cell, we calculated f karst as where i is one of the five categories of karstifiable rock and A overlay,i is the karstifiable area (km 2 ) in each 0.5°grid cell computed by overlaying the WOKAM polygons with grid cell polygons in the cylindrical equal-area projection.Share i is the fraction of each WOKAM polygon that is assumed to be covered by karstifiable rock.We arbitrarily set the Share to 0.4 for the discontinuous rock categories and to 0.9 for the continuous and mixed rock categories.f kmax is the maximum karstification fraction, which is accordingly set to 0.9.A cont is the WaterGAP continental area (km 2 ) per 0.5°grid cell.
The obtained f karst map (Figure 8d) reveals that China, Russia, the United States, and Canada have the largest karstifiable surface.The Middle Eastern countries (e.g., Iran, Egypt) and European countries (e.g., Montenegro, Bosnia and Herzegovina, Slovenia, France) have very high areal percentages of karst.In contrast, the fraction of karst in the southern hemisphere is low.Note that f karst indicates merely the surface-exposure of karstifiable rocks, so that some karst aquifers are not captured by WOKAM, including the Harmanköy-Beyyayla Karst system in Turkey (Aydin et al., 2013) and the Umm Er Radhuma aquifer in Saudi Arabia (Hoetzl, 1995).

Computation of GWR in Karst Areas
We developed the method for computing GWR in karst by comparing both GWR and runoff component R 3 as simulated by the standard version of WaterGAP 2.2e in cells with karst with ground-based estimates of karst GWR.Consistent with the knowledge that surface runoff is mostly negligible in karst areas (Hartmann et al., 2014), we found that a good fit to ground-based estimates is achieved if it is simply assumed that the total runoff from land, excluding urban runoff and soil overflow, recharges the aquifer in karst areas.This will be further discussed in Section 3.1.Accordingly, karst GWR is computed as R g karst is applied to the karst fraction of each grid cell, while GWR in the remaining area of the grid cell is computed with the method for non-karst (Equations 2 and 3).

Modifying the Computation of GWR Outside of Karst Areas
The global data sets of physiographic characteristics that are used to determine R gmax , f r , f t , f h , and f pg for each grid cell were produced in the early 2000s or even the last century (Döll & Fiedler, 2008).In this study, we implemented new global data sets of relief and soil hydraulic conductivity.We then optimized the relief-related factor f r , the soil parameters f t and R gmax (Figure 3) as well as the hydrogeology-related factor f h .Only the permafrost/ glacier-related f pg factor remains unchanged.Optimization occurred by maximizing the Nash-Sutcliffe efficiency (NSE, Equation 10) that results from comparing mean GWR as computed with a variety of plausible look-up table variants for computing the four parameters with the ground-based estimates of non-karst GWR in 422 0.5°grid cells.The optimized look-up tables are Tables A1-A3.

Water Resources Research
10.1029/2023WR036182 WAN ET AL.
The relief-related GWR factor is calculated based on the global EarthEnv 30 arc-second land surface elevation database (http://www.earthenv.org/topography),which was derived from the global 250 m GMTED2010 (Danielson & Gesch, 2011) and the near-global ∼90 m SRTM4.1dev(Jarvis et al., 2008).We adopted this database as it provides improved and spatially more resolved land surface elevation data compared to the 1 km GTOPO30 data set used to compute f r in the standard version of WaterGAP 2.2e (Döll & Fiedler, 2008); different from SRTM, it also includes high-resolution values north of 60°N.The median slope of the native DEMs within the 30 arc-second cells was downloaded.Eight slope classes from flat to steep are distinguished (Table A1) after checking against various slope classification schemes (e.g., Chang & Tsai, 1991;Fischer et al., 2000Fischer et al., , 2008;;Hammer et al., 1995).Overlaying the slope map onto the 0.5°grid-cell mask, the total area of each grouped slope class within each 0.5°cell was calculated, followed by the computation of the average relief r avg in each grid cell, with where A Cslope,j is the area of slope class j within the 0.5°cell, reflecting the distribution of slopes in one grid cell.
The relationship between r avg and f r shown in Table A1 reflects the expert knowledge that the fraction of runoff becoming GWR decreases with increasing slope and was determined by manually optimizing the fit of simulated GWR to the ground-based GWR estimates (jointly with optimizing heuristic relationship for soil-texture-related factors f t and R gmax in Table A2, and hydrogeological-related factor f h , see the next paragraph).The specific values for these factors are expert guesses and have been adjusted iteratively by comparing the resulting GWR.NSE was computed separately for humid and (semi)arid regions and the geometric mean of both was maximized.The f r for each cell is obtained by linear interpolation.The median grid cell values of r avg and f r are 2.55 (range 1-7.8) and 0.92 (range 0.22-1), respectively.
While up to this study, the soil-related factors f t and R gmax were based on information on soil texture (distinguishing the classes fine, medium and coarse) in the uppermost 30 cm of the soil from the FAO Digital Soil Map of the World and Derived Soil Properties (Döll & Fiedler, 2008;FAO, 1995), these two factors are now related to saturated hydraulic conductivity K sat of the soil (Figure 3b).The replacement is made because K sat is more directly related to the disposition of the soil to allow percolation of water to the groundwater than soil texture and a global-scale soil K sat data set is now available.Besides, soil K sat was found to be one of the most important predictors of natural groundwater recharge (Mohan et al., 2018;Racz et al., 2012).We used the global K sat data set (Montzka et al., 2017), with a spatial resolution of 0.25°and seven layers (0, 5, 15, 30, 60, 100, and 200 cm).This data set is an upscaling of the SoilGrids1km database (Hengl et al., 2014).For each 0.25°grid cell, we calculated the layer thickness-weighted harmonic mean K sat 0.25 based on the saturated hydraulic conductivity at z cm (0, 5, 15, 30, 60, 100) depth as 200 Overlaying this K sat 0.25 map onto the WaterGAP grid cell mask, areal average K sat for each 0.5°cell is obtained as the (arithmetic) mean K sat of four pixels.There remain around 8% of the WaterGAP land cells that do not overlap with any K sat 0.25 value (e.g., coastlines and islands); the global mean K sat (25.7 cm/d, log 10 K sat = 1.4) is assigned to these cells.The resulting K sat ranges between 1.1 and 1,168.1 cm/d with a mean of 49.9 cm/d.The 10th, 50th, and 90th percentiles of K sat are 8.0, 25.7, and 76.1 cm/d, respectively.As the global K sat shows a highly skewed distribution, the log scale of K sat , an alternative to better approximate a normal distribution, was used to interpolate soil-texture-related factors.
The groundwater fraction of total runoff in (semi)arid is expected to be lower than in humid areas for the same soil texture for several reasons (Döll & Fiedler, 2008), for example, due to surface crusting in the case of low vegetation cover and hydrophobic behavior of dry soil surfaces with organic material.Table A2 provides the optimized relationships between log 10 K sat and R gmax and f t , distinguishing (semi)arid and humid grid cells.The hydrogeology-related factor f h was increased as compared to the values of Döll and Fiedler (2008) used for WaterGAP 2.2e and is now computed using the look-up table Table A3.
In the case of all look-up tables, the grid cell values are computed by linear interpolation.The median values of R gmax , f t , and f h are 6.1 (range 3-7), 0.96 (range 0.5-1), and 0.9 (range 0.85-1), respectively.89% of 0.5°( semi)arid cells are constrained to the (semi)arid-coarse adjustment (Section 2.2), which has been modified as follows: if a (semi)arid grid cell has R gmax > 5 mm/d (indicating coarse texture), recharge occurs only when precipitation exceeds 12.5 mm/d.With the new values for f r , f t , and f h , the GWR factor f g was recalculated.The regional reduction of f g in Mississippi Embayment was removed because total water storage there as determined by GRACE shows an increasing trend while WaterGAP 2.2d showed a negative trend (Müller Schmied et al., 2021).

Comparison of Simulated GWR to Ground-Based Estimates
To optimize the GWR simulation approach, we aggregated the 1981-2010 daily R g (Equation 2) and R g karst (Equation 5) of the period 1981-2010 to represent the simulated long-term mean GWR (mm/yr).We decided to adopt a fixed 30-year period instead of a varying relevant period depending on the time stamp of observations for three reasons.(a) Only some ground-based estimates had information on the period of sampling, while most, for example, the estimates for Australia from Crosbie et al. (2010), did not; (b) The most widely used chloride mass balance reflects groundwater recharge over a long period before the measurement period; (c) The best groundbased GWR estimate for the 0.5°grid cell, which may be an aggregate of individual estimates from different time stamps or locations in the same grid cell (see Section 2.1), was used for comparison.

Computation of GWR per Grid Cell and the Contribution of Karst GWR to Total GWR in Regions
Daily GWR on land area per grid cell R g grid (mm/d) is calculated as the weighted sum of R g karst in the karst area and updated R g in the non-karst area, with where f land and f cont are the areal fraction of daily land area and continental area of the cell area; f karstland is the fraction of land area covered by karstifiable rocks, subject to the land area fraction.

Water Resources Research
10.1029/2023WR036182 WAN ET AL.
The regional contribution of karst GWR to total GWR is thus where j represents cells in the region of interest, A cell,j is the WaterGAP cell area within the selected region.

Performance Metrics
We used several performance metrics to compare simulated values of both GWR and streamflow Q to observations.The Nash-Sutcliffe efficiency (NSE) (Nash & Sutcliffe, 1970), which provides an integrated measure of model performance concerning mean values and variability and is computed as where μ obs is the mean of observations; sim(i) and obs(i) refer to the simulated and observed GWR respectively at the ith grid cell among n cells.It was used to evaluate the fit between simulated and observed long-term mean annual GWR, together with the coefficient of determination (R 2 ), the difference of the respective medians (dMD), and the percent bias (PBias), with The Kling-Gupta efficiency (KGE) together with its three components (Kling et al., 2012) is used to compare streamflow computed with the standard and the updated GWR approaches by WaterGAP to observations, with where CC is the correlation coefficient and RBias = μ sim μ obs (13) The three KGE components enable distinguishing the model performance regarding correlation, bias, and variability.KGE is applied for time series of monthly streamflow as well as for the statistical low flow Q 90 , the monthly streamflow that is exceeded in 9 out of 10 months.Streamflow observations are taken from Müller Schmied and Schiebener (2022).

Performance of Karst GWR Simulation
Only 8% of the 816 ground-based GWR estimates for 0.5°grid cells represent karst conditions; they are primarily located in southern Africa, northern China, and the Mediterranean (Figure 1).corresponding recharge ratio (long-term mean GWR to WaterGAP precipitation) can be as small as 2% and as high as 66%.Figures 4a and 5 show the comparison of long-term average GWR as simulated by the standard WaterGAP 2.2e model, which does not take into account karst and the updated GWR algorithm, to the groundbased estimates.The comparison reveals that the standard WaterGAP 2.2e model tends to strongly underestimate karst GWR.Underestimation is a common occurrence across almost all regions, with karst GWR in South Africa being underestimated by more than 90% (Figure 4a).Only 20% (13 out of 64) cells have GWR simulations that do not underestimate the observations by at least 30%.The median of the simulated karst GWR is, with 97 mm/yr, only half of the median of the observed GWR, and the regression curve of simulated versus observed GWR (green line in Figure 5a) shows an overall underestimation, too.This strong and expected underestimation of karst GWR by the standard WaterGAP 2.2e motivated the explicit consideration of karst in the GWR simulation by a new algorithm.Figure 4b displays the refined GWR with the updated GWR approach, which will be elaborated upon in the next two sections, specifically distinguishing karst and non-karst areas.
Figure 6 compares the 64 ground-based karst GWR estimates to three different sets of WaterGAP output: the standard GWR (R g in Equation 2), effective precipitation reaching the soil (P eff in Equation 1) and the major runoff component (R 3 in Equation 1).Different from standard GWR and P eff , WaterGAP 2.2e runoff R 3 achieves a rather unbiased estimation of ground-based estimates of GWR in karst, with a simulated mean of 270 mm/yr compared to the observed mean of 246 mm/yr, while the mean simulated WaterGAP 2.2e GWR for the 64 karst observations is only 101 mm/yr (Figure 6).However, the median value for R 3 is appreciably higher than the median of the ground-based estimates (277 vs. 197 mm/yr, Figures 4a and 5).If it is assumed that karst GWR is best simulated by R 3 , the NSE improves from a negative value for the standard WaterGAP 2.2e GWR to an NSE value of 0.55, which is not only due to the improved estimate of the mean but also an improved representation of the spatial variability of GWR (Figure 4); the coefficient of determination R 2 increases from 0.45 to 0.62 (Figure 6).In the 33 humid karst cells, GWR is slightly overestimated (+13.4% on average), while it is very slightly underestimated by 2.6% in the 31 (semi)arid cells.It is plausible to assume that runoff is a good estimator of GWR in karst as surface runoff is known to be very small in karst.Using the effective precipitation that reaches the soil as a predictor of GWR would lead to an overestimation of ground-based GWR by about a factor of 3, while the R 2 is between that of using standard GWR and R 3 as a predictor of karst GWR (Figure 6).With this new approach for estimating GWR in karst, the general underestimation of karst GWR in the Mediterranean is almost eliminated, and the significant underestimation in South Africa is decreased (Figure 4b).About half of the karst observations are simulated within 30%-50% of the (rather uncertain) ground-based estimates.Therefore, based on the comparison to ground-based GWR estimates, GWR is set to R 3 for the karst fraction f karst of each grid cell in WaterGAP (Equations 4 and 5).

Performance of Non-Karst GWR Simulation
Comparison of the standard WaterGAP 2.2e to the ground-based estimates in 752 non-karst grid cells shows that the model suffers from overestimation with small GWR values, approximately those below 10 mm/yr (Figure 5a) and underestimation with larger values (Figure S2c in Supporting Information S1).For model evaluation and development, we divided the grid cells with non-karst ground-based estimates into the 330 cells in Australia (AU) and the 422 cells outside of Australia (nonAU).WaterGAP tends to overestimate small GWR values, which are particularly distributed in Australia (Figures 4 and 5a).With an NSE of 0.50 and an R 2 of 0.52, model performance for the Australian non-karst cells is worse than for the non-Australia non-karst cells, which have an NSE of 0.56 (with an underestimation of the mean and median) and am R 2 of 0.64 (Figure 5a).Overestimation of GWR in Australia can be attributed to the specific native vegetation in Australia, which has developed deep root systems that can extract soil water at very low matric potentials and thus use rainfall more efficiently than vegetation elsewhere (O'Grady et al., 2008).GWR can be up to 100 times larger if native vegetation is replaced by annual crops (Crosbie et al., 2010), which explains that some of the high ground-based GWR estimates for Australia are  (Cleveland, 1979).
less or not at all overestimated by WaterGAP.These special characteristics of native vegetation, which lead to decreased runoff and GWR for the same land cover type as compared to outside of Australia, are not considered in WaterGAP, which uses MODIS land cover globally.Therefore, the Australian non-karst ground-based GWR estimates were not used to adjust the non-karst GWR simulation in WaterGAP.
Two-thirds of the non-Australian estimates are for (semi)arid cells (275) as GWR studies are more common in dry settings (Figures 5b and 5d).Considering the medians, GWR is underestimated in both climate zones by about 40%; lower observed values tend to be overestimated and higher observed values are underestimated in both zones (Figure 5b).In (semi)arid cells, underestimating runoff (points below the 1:1 line in Figure 5d and Figure S2b in Supporting Information S1) is partly responsible for the frequent underestimation of high GWR values, while in humid cells, the GWR underestimation appears to be more a result of an inadequate partitioning of runoff.Assuming that GWR is best simulated by R 3 , as for karst recharge, would lead to a strong overestimation of GWR, except for the high ground-based estimates in (semi)arid grid cells (Figure 5b and Figure S2b in Supporting Information S1).Despite the incommensurability and high uncertainty of ground-based estimates, it is worthwhile to try to improve the model fit to these estimates.
Adjusting GWR model parameters by optimizing the fit of simulated GWR to the non-karst ground-based GWR estimates outside of Australia leads to higher f g and R gmax values (Section 2.4 and Figures S3-S7 in Supporting Information S1) as well as to improved model performance.NSE increases from 0.56 for standard WaterGAP 2.2e to 0.65 for the updated GWR (Figure 7).The underestimation of the mean and the median of ground-based estimates is reduced, with the bias decreasing from 35% to 19%.R 2 has increased slightly from 0.64 to 0.67.Performance has improved in both humid and (semi)arid grid cells (Figure 7).The LOESS curve has become steeper, indicating less underestimation for high GWR values while maintaining performance for low GWR values.
For cells with GWR estimates in 5-250 mm/yr, there is no tendency to over-or underestimate GWR, with the regression curve close to the 1:1 line (Figure 7).The underestimates for GWR between 250 and 400 mm/yr are mainly located in Virginia, USA (Figure 4b for spatial comparison).Virginia is a humid subtropical region, and the 51 ground-based estimates derived from Sanford et al. (2012) show a total GWR of 37 km 3 /yr, which accounts for 78% of the simulated total runoff from land.The updated GWR is only 21 km 3 /yr but closer to the estimates than the standard simulation (16 km 3 /yr).For Australia, the updated GWR does not change too much, with the NSE dropping from 0.5 (Figure 5a) to 0.45 but R 2 slightly improving from 0.52 to 0.55 (not shown).Globally, the percentage of grid cells for which simulated GWR is within 0.7-1.5 times the ground-based estimate increased from 24% in the standard WaterGAP 2.2e to 31% in the updated model version (Figure 4).

Best Estimate of Global-Scale GWR and the Contribution of Karst GWR
Our current best estimate of mean annual GWR on the land area of each 0.5°grid cell (R g grid in Equation 8) during 1981-2010 is shown in Figure 8a.As compared to the standard WaterGAP 2.2e GWR, GWR is larger in the majority of grid cells, except in West Asia (Figure 8b).There are more grid cells with high GWR (>300 mm/ yr), while the distribution of low GWR (<2 mm/yr) is similar.The general increase is due to both the consideration of GWR and the updated non-karst GWR parameters.The annual variability in GWR is indicated by the coefficient of variation (CV), as shown in Figure 8c.An inverse relationship between GWR and CV is evident, where cooler, wetter regions typically exhibit lower variability (median CV of 25%), and drier regions have higher CVs, with standard deviations mostly exceeding 85% of the long-term means.
Based on WOKAM and our karst area parameterization (Section 2.3.1),we estimate that 11% (14.8 million km 2 ) of the global continental area (134.394 million km 2 excluding Greenland and Antarctica) is covered by exposed karstifiable rocks.Global renewable karst groundwater resources (excluding Greenland and Antarctica) are estimated to be 4,247 km 3 /yr, which amounts to 20% (C karst according to Equation 9) of the total renewable groundwater resources from diffuse groundwater recharge (Table 2).An earlier assessment using a simple budget equation estimated smaller karst groundwater resources of 3,165 km 3 /yr but a higher percentage of 26% (Stevanović, 2019).Compared to the standard WaterGAP 2.2e, with a total of 1,293 km 3 /yr of karst GWR taking into account the WOKAM-derived karstifiable area to compute karst GWR, the total GWR in karst areas has almost tripled, resulting in a doubled contribution to total GWR (Table 2).
Karst GWR accounts for nearly one-third of GWR in North/Central America and Asia but for only 6% in South America (Table 2).In (semi)arid cells, the contribution is, with 30%, on average higher than in humid cells, with 20%, even though humid cells have a larger extent of karst areas (Figures 1 and 8d).This disparity is due to the lower f t values (Table A2) and the (semi)arid-coarse adjustment (Section 2.4) applied to the majority of (semi)arid cells and; the latter adjustment restricts recharge in the non-karst fraction when precipitation is below 12.5 mm/d, whereas it is assumed that this restriction does not apply in the karst fraction.Nearly one-third of global cells in over 140 countries have non-zero f karst (Figure 8d), with 65% of these cells being humid.We categorized the 0.5°c ells into three classes: those cells with f karst less than 0.05, between 0.05 and 0.5, and larger than 0.5.In these three classes, 0.8% of the precipitation on land (and 1% of total runoff), 10% of the precipitation (24% of the runoff), and 30% of the precipitation (76% of the runoff) recharge the respective karst aquifers.GWR as a fraction of total runoff from land (GWRF) is shown in Figure 8e.The cells with high cell GWRF (>0.7) generally coincide with high f karst (Figure 8d).For all land areas excluding Antarctica and Greenland, the long-term mean GWR, that is, the renewable groundwater resources, is estimated to be about 21,000 km 3 /yr, which is 56% higher than the estimate of the standard WaterGAP 2.2e (Table 2).According to our updated estimation, GWR accounts about half the total renewable water resources (from diffuse recharge) of the Earth, and not only for about one-third as in the standard WaterGAP 2.2e.39% of this increment is due to the explicit consideration of karst, while the rest is the result of the modified GWR approach for non-karst areas, in particular the increased groundwater recharge factor (Figure S7 in Supporting Information S1).Continental long-term average GWR values range between 446 km 3 /yr for Australia and Oceania to 6,477 km 3 /yr for Asia, and the renewable groundwater resources to renewable water resources fraction is, with 36%, lowest in Australia, and, with 61%, highest in Europe (Table 2).Considering the gridded human population in 2020 (Rose et al., 2021), the per-capita renewable groundwater resources of continents are generally negatively correlated to the population density (Table 2).Coefficient of variation of (updated) annual GWR (c), which is the standard deviation of annual GWR to the mean.Karstifiable fraction f karst of the continental area of the 0.5°grid cell in % (d).Updated GWR as a fraction of total runoff from land (e).

Validation of Simulated GWR Using Streamflow Observations
GWR affects the dynamics of streamflow as groundwater discharge to surface water bodies is delayed as compared to surface runoff.Different from ground-based GWR estimates, which are (mainly) point estimates, streamflow observations reflect the aggregated processes in the upstream drainage basin, and thus the GWR within the whole basin.With increased GWR, it is expected that the variability of streamflow decreases as increased groundwater discharge to streams during dry periods leads to an increase in the low flows.Therefore, the validity of the new algorithms for GWR computation in WaterGAP is tested by comparing the fit of simulated and observed time series of monthly streamflow for both the standard WaterGAP variant and the variant with the updated GWR algorithms.Only the standard model variant was calibrated against observed mean annual streamflow, while the calibration parameters of the standard variant were also used for the variants with updated GWR.With the updated model and the consequently strongly increased GWR, the already underestimated temporal variability of streamflow at 1,509 gauging stations is, on average, underestimated even more (RVar in Figure 9a), while the correlation is only slightly decreased (CC in Figure 9a).Evaluating the fit of monthly low flow Q 90 for a subset of 1,357 gauging stations with more than 15 years of observations and a drainage area of more than 20,000 km 2 shows that Q 90 is overestimated more strongly with the updated GWR than with the standard GWR (Figure 9b).As the impact of GWR on streamflow is affected by the choice of the groundwater discharge coefficient, with groundwater discharge to the surface water bodies being the product of groundwater storage and discharge coefficient, a higher coefficient delays the response of groundwater to GWR less strongly.If the globally constant discharge coefficient is increased from 0.01/d to 0.1/d, the updated GWR leads to an unbiased streamflow variability and a slightly improved correlation as compared to the standard WaterGAP 2.2e variant (Figure 9a).The fit of the simulated Q 90 values improves as compared to the standard variant, with KGE increasing from 0.81 to 0.93 (Figure 9b).

Sensitivity of Estimated GWR to Climate Data Set
The fifth generation of the ECMWF reanalysis (ERA5) and its predecessors are one of the most popular atmospheric reanalyses by the climate community (Dee et al., 2011), for example, for the evaluation of CMIP6 (Harvey et al., 2020).The W5E5 data set used in this study is a bias-adjusted version of ERA5, where, among other variables such as shortwave radiation, daily precipitation was adjusted using the observed monthly precipitation totals and the observed number of wet days per month (Cucchi et al., 2020;Lange et al., 2021).We applied updated GWR algorithms for both karst and non-karst areas to the runoff from land of a version of WaterGAP 2.2e that was driven by GSWP3-ERA5 and calibrated against observed streamflow in the same way as the version Note.For the updated algorithm only, the total runoff from land (B) as well as the GWR as a fraction of total runoff from land (E/B), and the per capita GWR (E/A) is listed, together with the population (A). a Excluding Antarctica and Greenland.
b Population data based on gridded LandScan database for the year 2020 (Rose et al., 2021).c C karst according to Equations 8 and 9, where karst recharge in the standard computation is the standard GWR while in the updated computation is runoff R 3 .d Per capita GWR per day.e Eurasia is subdivided into Europe and Asia along the Ural, and Turkey is assigned to Asia.

Water Resources Research
10.1029/2023WR036182 WAN ET AL.
driven by GSWP3-W5E5 (Section 2.2). Figure 10 shows the difference between GWR and precipitation, which reflects the sensitivity of the GWR algorithm to precipitation.
The global renewable groundwater resources calculated with WaterGAP 2.2e under GSWP3-ERA5 is 23,000 km 3 /yr, which is higher than that under GSWP3-W5E5 (Table 2).While the spatial distributions of GWR under both data sets are very similar (not shown), the differences in GWR between the two climate data sets correspond closely to the variations in precipitation (Figures 10a-10c).ERA5 data indicates a 7% higher global average precipitation, especially in China and western South America; in almost one-sixth of the 0.5°grid cells, mean precipitation differs by more than 50% (Figure 10c).The higher global precipitation estimate, combined with adjustments made to the runoff coefficient by model calibration among model runs, leads to a similar increase of simulated GWR.
When compared to the 816 ground-based estimates, GSWP3-ERA5-driven GWR performs better in terms of capturing the mean and median than GSWP3-W5E5; however, it exhibits a lower correlation and NSE value (Figure 10d).The ERA5 precipitation in these 816 cells is on average 4% lower than the W5E5 precipitation,

Sensitivity of Simulated GWR to Karst Area Parametrization
WOKAM does not provide exact information on karst (or rather karstifiable area) but just rough classes; for example, a map unit labeled "discontinuous karst" means that 15%-65% of the map unit area is karstifiable, and we arbitrarily set the area share value for discontinuous karst to 0.4 (Section 2.3.1).To evaluate the impact of karst area uncertainty on global karst area and groundwater recharge, we carried out simulations with the parameter area share in Equation 4 varying from 0.15 to 0.65, at intervals of 0.05, for discontinuous rock categories, and from 0.7 to 1.0 for continuous and mixed rock categories.The karst recharge and its contribution to global total GWR are derived from Equations 8 and 9.
According to this sensitivity analysis, global karst(ifiable) area as mapped in WOKAM may cover 8%-13% of the global land area (excluding Greenland and Antarctica), compared to 11% for the two baseline values applied in this study (Figure 11a).The contribution of karst recharge to global GWR may vary between 14.7% and 23.5%, with the baseline value being 20% (Figure 11c).Total global GWR could range from 20,450 to 21,670 km 3 /yr, that is, by 3.4% to +2.3% of the baseline value (Figure 11b).However, this variance is much smaller than the impact of climate input.Using the GSWP3-ERA5 data set could increase global GWR to 23,000 km 3 /yr, that is, by +8.5% of the baseline value driven by GSWP3-W5E5 (Section 4.1).

Comparison to Previous Studies
Large-scale studies on karst GWR are rare.Hartmann et al. (2015) presented the karst recharge model VarKarst-R for karst regions in Europe and the Mediterranean, distinguishing four karst landscapes with different parameter sets.They compared the results of both their model and a WaterGAP version from the year 2008 to 36 groundbased estimates of karst GWR.They found that different from VarKarst-R, this WaterGAP version, which did not take into account karst, significantly underestimated the ground-based estimates (Figure 12a).A comparison of both standard and updated GWR estimates of WaterGAP 2.2e to the same ground-based estimates shows that the standard WaterGAP 2.2e still significantly underestimates the ground-based estimates (blue points in Figure 12b).WaterGAP 2.2e with the updated GWR algorithm no longer tends to an underestimation, and, with an NSE of +0.27, even outperforms VarKarst-R, which results in a negative NSE of 0.33.This demonstrates the effectiveness of the updated GWR algorithm, at least for Europe and the Mediterranean.
In addition, we compared the updated GWR to the data used in Berghuijs et al. (2022), global database of GWR at 5207 (most likely non-karst) locations compiled by Moeck et al. (2020) and found that the correlation has improved due to the updated GWR approach from 0.55 to 0.66 and the underestimation of the mean has decreased (Figure S8 in Supporting Information S1).The underestimation in the wet landscapes of Africa found in West et al. ( 2023) is strongly reduced (Figure 4), where the median discrepancies in annual recharge for these 32 wettropical points with the standard WaterGAP 2.2e and the updated version are 28.3 and 11.7 mm/yr, respectively.
For a nation-scale assessment, we compared the updated GWR to the national estimates available for the seven European countries (UK, Ireland, Finland, France, Spain, Netherlands, and Denmark).These estimates were  obtained through different approaches and at varying resolutions (Martinsen et al., 2022, their Table 1).We aggregated them to the spatial resolution of 0.5°and the resulting gridded mean annual GWR values are mapped in Figure 13a.The average GWR for the seven countries is 200 mm/yr (344 km 3 /yr) according to the national estimates and 260 mm/yr (446 km 3 /yr) according to the updated WaterGAP.Except for Spain, the updated WaterGAP results in higher GWR estimates than the national sources.The simulated GWR for countries other than Ireland and Finland show close agreement, with bias ranging between 18% for Spain and 34% for Denmark (Figure 13c).In Ireland, despite the high precipitation and a high presence of karst (Figure 8d), the national GWR estimate is, with 167 mm/yr, surprisingly low, and much lower than our simulated value of 598 mm/yr.This low recharge might be attributed to the threshold approach they used, where aquifers are assumed to have low recharge acceptance capacity (Williams, 2011).In Finland, GWR was estimated using the same threshold approach, which is to represent the presence of low-permeability bedrock at a very shallow depth below the land surface, which results in fast surface and subsurface runoff that flows directly to the surface water system.Although the low permeability of the aquifers in Finland relative to other world regions is captured by WaterGAP (Figure S6a in Supporting Information S1), we overestimate recharge in Finland by 80% because the K sat values according to Montzka et al. (2017) are high for Finland (Figure S4 in Supporting Information S1).Overall, when evaluating the grid-cell GWR, the performance is acceptable, with the regression curve close to the 1:1 line and an R 2 of 0.45.

Conclusions
The approach for computing diffuse groundwater recharge (GWR) in the global hydrological model WaterGAP has been improved by optimizing, for the first time, simulated GWR against a large number of ground-based GWR estimates, which were aggregated to best ground-based estimates in the 0.5°grid cells of WaterGAP.Also for the first time, GWR in karst areas has been taken into account explicitly in a global hydrological model, utilizing information on karstifiable areas from the WOKAM map.Based on ground-based estimates of karst GWR in 64 0.5°grid cells, we found that total runoff from land (except overflow and urban runoff) is best assumed to recharge the groundwater in karst areas.This results in an increase of NSE from 0.14 to 0.55, which is not only due to an increase in the mean but also an improved representation of the spatial heterogeneity.In addition, the fit of simulated GWR to ground-based estimates is similar to the fit achieved by the European-scale karst recharge model of Hartmann et al. (2015).
The fit to ground-based estimates of non-karst GWR in grid cells was improved by using up-to-date relief and soil information and by adjusting parameters and algorithms based on ground-based estimates for 422 grid cells outside of Australia.We found that assuming a globally homogenous groundwater recharge factor f g could lead to even higher NSE for the 814 0.5°grid cells with ground-based GWR estimates.We attribute this and any discrepancies between simulated GWR and ground-based estimates mainly to (a) the spatial incommensurability of ground-based estimates and the spatial extent of the WaterGAP grid cells (∼2,500 km 2 at the Equator) as plotscale measurements within such an extent can hit a large variety of GWR values and (b) deviation of the simulated total runoff from reality.In addition, the time period of ground-based observations was not considered in the study as this information is very often not available, which might lead to large discrepancies due to interannual climate variability.We are confident that physical properties such as karst features, relief, soil hydraulic properties, hydrogeology, and the occurrence of permafrost and glacier do impact GWR in addition to precipitation, AET, and total runoff (e.g., Moeck et al., 2020;Mohan et al., 2018;Scanlon et al., 2010;West et al., 2023) and that it is, therefore, appropriate to take these properties into account when quantifying GWR.The application of two alternative global climate data sets leads to similarly good model performance with respect to ground-based GWR estimates.A more detailed global sensitivity analysis of the parameterizations of physical properties will be conducted in the future to enhance our understanding of the model results.
With the updated GWR approach, GWR estimates are considerably higher than with the standard approach almost everywhere.For all land areas excluding Antarctica and Greenland, the renewable groundwater resources from GWR are estimated to be 21,000-23,000 km 3 /yr, 20% of which is generated in karst areas.The plausibility of the comparably large simulated GWR values is difficult to judge.While the simulated GWR values still tend to underestimate ground-based estimates unless mean GWR is more than 200 mm/yr, independently calculated GWR estimates for seven European countries are smaller than our model values, except for Spain.This points to the heterogeneity of ground-based and national modeling-based values.Assessing the updated GWR approach by comparing streamflow as simulated by WaterGAP 2.2e with the standard and the updated approach to observations showed that GWR cannot be validated by streamflow observations only, as simulated groundwater discharge and thus streamflow dynamics depend not only on GWR but also the assumed groundwater discharge coefficient.However, with the new GWR and a globally increased groundwater discharge coefficient, an unbiased temporal variability of monthly streamflow could be achieved for 1,509 streamflow gauging stations, and the computation of low flows was improved as compared to the standard WaterGAP 2.2e.
With the presented improved but still uncertain global-scale estimation of GWR, renewable groundwater resources account for more than half of the total renewable water resources, while with the standard approach, it accounted for only one-third.Given that groundwater is better protected from pollution and provides storage for the temporally highly variable runoff, the increase in GWR estimates is good news.We recommend that developers of global hydrological models utilize the WOKAM map of karstifiable rocks as well as ground-based estimates of groundwater recharge to improve the simulation of diffuse groundwater recharge in their models.

Appendix A
Look-up tables for factors in the groundwater recharge computation (Tables A1-A3).Note.The slope is expressed as a percentage, not as an angle because laboratory experiments that analyze the relationship between infiltration and slope gradients use slope mostly in the form of % due to the ease of measurement.

Figure 1 .
Figure 1.Spatial distribution of the ground-based GWR estimates for 0.5°grid cells in the compiled GWR data set.The locations of 64 grid-level estimates of karst GWR and 752 estimates for non-karst GWR are shown on a background of humid and (semi)arid regions as defined in WaterGAP.

Figure 2 .
Figure 2. Schematic of diffuse GWR simulation driven solely by precipitation in WaterGAP 2.2e.Boxes represent water storage compartments; arrows represent water flow, and (…) indicates other lateral water balance processes that do not affect GWR simulation.GWR simulation in model version 2.2e is the same as in 2.2d, which was described in detail by Müller Schmied et al. (2021).

Figure 3 .
Figure 3. Workflow of mapping the 0.5°relief-related parameter f r (a) and of the soil texture-related parameters f t and R gmax (b).

Figure 4 .
Figure 4. Ratio of simulated GWR to ground-based GWR estimates in karst cells (triangles) and non-karst cells (dots).Simulated GWR in (a) is the output of the standard WaterGAP 2.2e, while simulated GWR in (b) is the updated GWR.Red and yellow denote a GWR underestimation by WaterGAP and blue overestimation.

Figure 5 .
Figure 5.Comparison of ground-based GWR estimates to simulated values of GWR (a, b) and runoff R 3 (c, d) per grid cell as computed by the standard WaterGAP 2.2e, plotted in the log-log scale.The results are grouped according to estimates in non-karst areas in Australia (AU) and outside of Australia (nonAU) and in karst areas (a, c), and, only for non-karst cells outside Australia, grouped according to (semi)arid and humid conditions (b, d).The dashed lines are the LOESS (nonparametric locally weighted least squares regression) curves(Cleveland, 1979).

Figure 6 .
Figure 6.Comparison of ground-based karst GWR to R g , effective precipitation P eff , and runoff R 3 as simulated by WaterGAP 2.2e.μ obs and μ sim are the mean of ground-based estimates and simulates for the respective variables.

Figure 7 .
Figure 7.Comparison of ground-based non-karst GWR outside of Australia with standard WaterGAP 2.2e GWR and GWR as updated in this study.The black and red lines are the LOESS curves for the standard and updated GWR, respectively.The points are grouped according to (semi)arid and humid conditions.

Figure 8 .
Figure 8. Updated estimates of long-term average diffuse groundwater recharge GWR per grid cell for the period 1981-2010 (a).Change of GWR in mm/yr compared to standard GWR (updated minus standard) (b).Coefficient of variation of (updated) annual GWR (c), which is the standard deviation of annual GWR to the mean.Karstifiable fraction f karst of the continental area of the 0.5°grid cell in % (d).Updated GWR as a fraction of total runoff from land (e).

Figure 9 .
Figure 9. Performance of the simulated time series of monthly streamflow (a) and the simulated statistical low flow Q 90 (b) of three WaterGAP model variants, with the standard GWR algorithm (standard GWR, white symbols), the updated GWR algorithm (updated GWR, gray symbols), and with updated GWR algorithm and a tenfold higher groundwater discharge coefficient in all grid cells (updated GWR and discharge coeff., green symbols).(a) KGE and its three components correlation, bias, and variability (Equations 12-14) for the monthly streamflow time series at the 1,509 WaterGAP 2.2e calibration stations, (b) comparison of simulated and observed monthly Q 90 at the 1,357 calibration stations with more than 15 years of observations and a drainage area of more than 20,000 km 2 .

Figure 10 .
Figure 10.Difference of mean annual GWR computed with the updated GWR algorithm in mm (a) and in % (b) between WaterGAP 2.2e driven by the GSWP3-ERA5 climate forcing and the GSWP3-W5E5 climate forcing.The relative difference (%) of annual mean precipitation between the two climate forcings (c).Comparison of ground-based estimates (both karst and non-karst estimates including Australia) and updated GWR computed with climate forcings (d).In the case of positive values in (a)-(c), the GSWP3-ERA5 climate forcing leads to higher values.

Figure 11 .
Figure 11.A comparison of simulated global total karst area (a), global total GWR (b), and karst recharge contribution to total GWR (c) across varying karst area shares for karstifiable rock categories.The blue dots represent the results under the baseline area share values.

Figure 12 .
Figure 12.Comparison between 36 karst recharge estimates used in Hartmann et al. (2015) and their model result (VarKarst-R) and three versions of WaterGAP output.

Figure 13 .
Figure 13.Mean annual groundwater recharge at 0.5°spatial resolution according to national estimates of seven European countries, in mm/yr (Martinsen et al., 2022) (a).GWR as computed by WaterGAP 2.2e with the updated GWR approach (b).Comparison of gridded GWR from national estimates and WaterGAP (sim) (c).Orange colors in c indicate more data points (grid cells).The mean GWR according to national estimates (obs) and WaterGAP (sim) per country are listed in the table in the lower corner of panel (c).

Table 1
Summary of Data Sets and Additional Karst Studies Used for Compiling Global-Scale Ground-Based GWR Estimates

Table 2
Continental and Global Sums of Total (Diffuse) GWR (C and E) and the Fractions Stemming From Karst Recharge (D and F) as Computed by WaterGAP Using Both the Standard and the Updated GWR Algorithms for thePeriod 1981Period  -2010

Table A1 Slope
Classes and the Relief-Related Groundwater Recharge Factor f r

Table A2
Relationship Between the Saturated Hydraulic Conductivity of the Soil K sat (cm/d) and the Maximum Groundwater Recharge Rate R gmax and the Soil Texture-Related Recharge Factor f t Average annual temperature >15°C and average annual precipitation >1,000 mm in the period1961-1990. a