Groundwater flow across spatial scales: importance for climate modeling

Current regional and global climate models generally do not represent groundwater flow between grid cells as a component of the water budget. We estimate the magnitude of between-cell groundwater flow as a function of grid cell size by aggregating results from a numerical model of equilibrium groundwater flow run and validated globally. We find that over a broad range of cell sizes spanning that of state-of-the-art regional and global climate models, mean between-cell groundwater flow magnitudes scale with the reciprocal of grid cell length. We also derive this scaling a priori from a simple statistical model of a flow network. We offer operational definitions of ‘significant’ groundwater flow contributions to the grid cell water budget in both relative and absolute terms (between-cell flow magnitude exceeding 10% of local recharge or 10 mm y−1, respectively). Groundwater flow is a significant part of the water budget, as measured by a combined test requiring both relative and absolute significance, over 42% of the land area at 0.1° grid cell size (typical of regional and mesoscale models), decreasing to 1.5% at 1° (typical of global models). Based on these findings, we suggest that between-cell groundwater flow should be represented in regional and mesoscale climate models to ensure realistic water budgets, but will have small effects on water exchanges in current global models. As well, parameterization of subgrid moisture heterogeneity should include the effects of within-cell groundwater flow.


Introduction
Groundwater comprises Earth's main store of liquid freshwater. As such, groundwater has emerged as a critical source of water for irrigation and other uses. Groundwater quantities and flows are usually difficult to observe [1]. However, progress has recently been made in quantifying steady-state global groundwater flow patterns [2] and human withdrawals [3][4][5]. Local and regional studies have also emphasized the role of groundwater in sustaining ecosystems ranging from the Amazon rainforest [6,7] to desert oases [8,9], a role vulner-Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. able to human water diversions and to anthropogenic climate change [10][11][12][13].
Groundwater flow patterns are controlled by topography, geology, climate, and vegetation [14,15]. Groundwater sustains surface water bodies during drought periods, and may flow across topographically defined watershed boundaries, so that the nominal upstream drainage area need not represent the actual contributing area for a surface water body [16,17]. While groundwater flows and their coupling with the soil and atmosphere have been modeled at the scale of individual watersheds [18], it is not yet clear how to represent the effects of these flows on budgets at larger scales, for example at the grid resolution of global climate models [19]. Isotope and ion tracers have been used to infer that groundwater flow between Costa Rican watersheds with length scales ∼1 km can account for over half of streamflow [20]. Based on a water balance approach which compared measured streamflow to estimated basin groundwater recharge, it has been argued that many USA watersheds with areas typically 10 2 -10 5 km 2 export or import appreciable fractions (∼10-100%) of their recharge [21]. On the other hand, modeling of groundwater flow over the Danube basin in central Europe at 5 km grid resolution found that flow between grid cells was quantitatively insignificant in comparison to recharge, with typical magnitudes well under 1 mm y −1 [22].
Here, we estimate between-cell flows at different grid aggregation levels, and map their importance globally, using results from an equilibrium groundwater model run with topography resolved to 30 (0.0083 • or ≈1 km) [2]. We also estimate the fraction of grid cells where groundwater flow is a significant component of the water budget at different grid resolutions. We compare our results on the scaling of between-cell groundwater flows with grid spacing with those we obtain from a conceptual statistical model of groundwater flow in an idealized landscape. Although the significance of lateral groundwater flow depends on climate, topography and geology-e.g., more recharge, steeper terrain, and thicker aquifers favor greater lateral groundwater convergence [21]here we focus on the effect of spatial scale, in the hope of deriving a simple scaling relationship that can directly inform the development of the land component of climate and Earth system models.

Global numerical model simulations
The groundwater budget for any land area, such as a model grid cell, may be written aṡ It includes the rate of change in water contentV , recharge R, discharge D, and net lateral flow from neighboring areas L, all expressed as water amount per unit area per unit time (e.g. mm y −1 ). L may be positive or negative. At steady state,V = 0. Artificial extraction of groundwater and plant uptake of shallow groundwater may be considered as part of the discharge term.
We simulated steady-state groundwater flow globally (Greenland and Antarctica excepted) with modern-day sea level as the hydraulic head boundary condition. Details of model formulation and validation are described in [2], and here we only outline the key concepts. The steady-state model solves the mass balance equation within each grid cell and calculates inter-cell flow by Darcy's law. Groundwater recharge (R) is from a global land surface model simulation by [23] validated with >2000 streamflow observations worldwide and shown to best reproduce water table and wetland observations [2]. This recharge is redistributed laterally (L) from high to low water table elevations by Darcy's law. Thus, in upland cells the recharge (R) balances lateral divergence (negative L), and in valley and coastal cells convergence (positive L) causes the water table to rise to the land surface, which triggers a modeled groundwater sink as surface drainage or evaporation (D). Permeability was assumed to decline exponentially with depth, with the surface values depending on mapped soil type and the rate of decrease with depth depending on land slope (flat valleys accumulate deep sediment so are expected to have greater permeability at given depth). We performed the simulation continent by continent at 30 resolution (≈1 km at mid-latitudes). Starting the initial water table at the land surface, we solved the equations iteratively until the deviation from steady state,V , is below 1 mm y −1 everywhere. Cell-to-cell groundwater flow was recorded at the end of the simulation, and the net lateral flux L computed for each cell.

Calculation of flow scaling
We aggregated L from the original 30 cells into nonoverlapping 2 × 2, 3 × 3, . . . windows. Global area-weighted mean |L| were computed for aggregation levels ranging from the original 30 (1 × 1 windows) to 10 • (1200 × 1200 window size), spanning the range of grid resolutions in regional and global climate models. Only aggregated cells with over 50% land area were included in the global means. We highlight the results for 0.1 • (≈10 km) resolution, representing a typical grid for current to near-future regional and mesoscale climate models, and for 1 • (≈100 km) resolution, representing a typical grid for global climate models [24][25][26]. We show maps of L , R, and modeled water table depth at the 1 • resolution to give a sense of the geographic distribution of positive and negative L.
We also quantified the fraction of land area for which lateral flow L could be considered significant in magnitude at different grid resolutions. We defined a relative significance threshold as |L| > 0.1R, implying that lateral flow makes a substantial contribution to the water budget relative to the recharge term. However, in arid areas where R is small (e.g. ∼1 mm y −1 ), this criterion may be met with small absolute values of L. We also defined an absolute significance threshold as |L| > 10 mm y −1 , a value which approximates 0.1 of the globally typical recharge (mean R based on [23] was 107 mm y −1 ). In very moist areas, this may be small relative to local recharge. Therefore, our combined criterion for significance required both |L| > 0.1R and |L| > 10 mm y −1 .
Note that unlike the lateral flow magnitude |L|, which as we shall see is strongly scale dependent, the expectations of the recharge R and discharge D terms in equation (1) are not inherently scale dependent, though they may vary as a function of, for example, climate and geology. Intuitively, this is because both these quantities are defined to be nonnegative; by contrast, L may be of either sign, and as we average across many local groundwater exporting and importing areas with respectively negative and positive L, it becomes a smaller term in the water budget relative to R and D. Thus, while average |L| is expected to decrease with increasing grid spatial scale, |R| and |D| are not expected to strongly vary with grid scaling. Over a region large enough that L may be neglected, equation (1) indicates that steady-state mean groundwater discharge D is essentially equal to mean groundwater recharge R.

A statistical model for groundwater flow scaling
To gain more insight into the dependence of groundwater flow on grid spacing, we developed a simple statistical model whose predictions as to the scaling of mean |L| can be compared with our detailed numerical simulations. This model is based on an idealized steady-state landscape comprised of many pixels, with the groundwater budget in each pixel obeying equation (1). The recharge rate R is taken to be constant within the landscape, while discharge is confined to an area fraction f and uniformly distributed across pixels within this fraction. For real landscapes, typical discharge area fractions f , which may comprise concentrated or diffuse discharge features such as springs, seeps, and stream channels, have been estimated to range from ∼0.01 for quite arid areas to ∼0.5 for wetland-dominated areas with high recharge rate, low relief, and low permeability [27,8,28,29].
Given the assumption of steady state, groundwater must travel from source areas (1 − f of the landscape, recharge only) to sink areas (with both recharge and discharge). L at the source pixels will be equal to −R, while we set L at the sink pixels to 1− f f R to maintain mass balance at the landscape level: groundwater flows across the landscape boundaries are assumed negligible, so L averaged across the landscape must tend to zero. That is, If the grid scale is fine enough that it fully separates source from sink areas, the average lateral flow magnitude (i.e. the mean absolute value of L) will be an area-weighted average of that from source and sink pixels: where · denotes an expectation taken across the landscape grid cells and the subscript 0 denotes a resolution fine enough to resolve source and sink pixels. Assuming f < 0.5, this means that at fine enough grid resolution the average lateral flow magnitude will be somewhat larger than the recharge magnitude. However, at larger grid scales, source and sink areas will be found within the same grid cells, which will result in smaller |L| at the grid cell level. If the sink area fractions within each grid cell are represented by φ (which has an average value of f but will differ between grid cells), then the lateral flow magnitude averaged over these larger grid cells will be This lateral flow magnitude per unit area is again simply an area-weighted average of the lateral flow in source areas (−R) and that in sink areas ( 1− f f R). Suppose the typical length scale of a discharge area is s and the grid cell length scale is S. Regard sink areas as pixels randomly and independently distributed with respect to grid cell boundaries, with each grid cell containing exactly N = (S/s) 2 total pixels, each of equal area. Within each grid cell, let n designate the number of sink pixels (leaving N − n source area pixels). n would then follow a binomial distribution with parameters N and f . The mean absolute deviation of this distribution, MAD = |n − f N | , can be shown [30] to be where ν is the least integral greater than f N and P ν = ( N ν ) f ν (1 − f ) n−ν is the probability of the outcome ν under this binomial distribution. Since φ = n/N , the expected value of | φ f − 1| would be equal to MAD f N . Substituting into equation (3), Note that this expression reduces to equation (2) for fine grids (S = s or N = 1, in which case ν = 1 and P ν = f ). In the limit of large N , ν → f N and Substituting these values into equation (5), we get for large grid spacings Thus, based on this simple statistical model we expect mean lateral flow |L| to scale as the reciprocal of grid length scale S.
We applied equation (5) over the globe, with R being the global mean recharge from [23] (the spatially varying version is used to drive the numerical model discussed in the previous subsections). The numerical model also provided an estimate of f , the global discharge area fraction, while s is taken to be the model resolution of 30 . We compared this statistical scaling model with the results of aggregating our numerical simulation of groundwater flow to estimate |L| at different grid spacings.

Flow scaling
Global mean lateral flow magnitude |L| drops steadily with increasing grid spacing, going from 167 mm y −1 (156% of mean recharge) at the original resolution of 30 to 24 mm y −1 (22% of mean recharge) at a resolution of 0.1 • and 2.2 mm y −1 (2.0% of mean recharge) at a resolution of 1 • (figure 1(a)).
The area fraction of grid cells with significant lateral flow also drops with aggregation to coarser grid resolutions, as expected ( figure 1(b)). The fraction of grid cells meeting our relative criterion decreases from 99% at the original resolution to 78% at 0.1 • and 26% at 1 • . For the absolute criterion the respective figures are lower, 79%, 47%, and 3.2%, and for the combined criterion lower yet, at 79%, 42%, and 1.5%.

Mapping of lateral flow
At the 1 • resolution, lateral flow can be visualized on a global map ( figure 2(b)). L has occasional relatively large positive values in the neighborhood of 10 mm y −1 (indicating groundwater inflow) in grid cells across a variety of geographical and climatic settings that receive groundwater convergence, from desert areas (for example in Australia and Arabia) to moist tropical and temperate regions (e.g. in Amazon Basin and Poland). Particularly in arid areas where recharge is low (figure 2(a)), groundwater inflow of this magnitude could be important in sustaining relatively shallow water tables (figure 2(c)) and groundwater-fed surface water features, though these will typically occupy only a fraction of a 1 • grid cell. Some substantial negative values of L (indicating outflow) occur in steep coastal areas, such as in Norway, where there is flow directly into the ocean, and to some extent in highland blocks (Appalachians, central New Guinea) where there is net flow toward valleys. Consistent with Darcy's law, larger values of |L| in figure 2(b) tend to occur in areas where the water table height (surface elevation minus the water table depth shown in figure 2(c)) shows large changes between cells.
At 0.1 • resolution (figure 3), lateral flows are larger by an order of magnitude than at 1 • , as expected, and widespread flows of order 50 mm y −1 occur, mostly within 1 • grid cells, in moist areas with moderate relief such as Europe and southeast Asia. At this resolution more oases can also be seen in arid areas such as the Sahara desert.

Comparison with expectations from statistical model
We found that the mean lateral flow magnitude at our original grid scale to be 1.56 times the global mean recharge R. This in reasonable agreement with our expectation from the simple statistical model presented above given the global discharge area fraction f = 0.144 found in our model: equation (2) would give |L| = 1.71R at the model grid scale, some 10% higher than observed. The modest discrepancy may be attributed to the tendency of f to be positively correlated with R (more extensive discharge where climate is wetter), so that the recharge-weighted discharge area fraction f that should be used in equation (2) would be larger than the global mean value.
Global mean lateral flow magnitude |L| is very close to proportional to the reciprocal of the grid spacing (as expected from equation (6)) over a broad range of grid scales, between about 0.02 • and 2 • (figure 1(a)). Close inspection of figure 1(a) shows that over this range, |L| exceeds by about 30% the quantitative predictions made with equation (5) using the global mean f and R and with s taken equal to our model grid scale of 30 , perhaps because the effective mean discharge area size in our numerical simulations was somewhat larger than one 30 cell.
|L| decreases more slowly than the reciprocal of the grid spacing close to the original resolution (grid resolutions under 0.02 • ). This may be because this length scale is still close to the typical discharge area size in our groundwater model.
|L| also decreases more slowly than the expected dependence on grid scale at resolutions coarser than 2 • , by which point the average |L| has dropped below 1 mm y −1 ( figure 1(a)). This may reflect (1) long-distance flow between climate zones with differing recharge that is not captured in the conceptual framework; (2) flow directly into the sea from coastal grid cells, which does not cancel as neighboring land grid cells are aggregated; or (3) numerical error in our groundwater model, which was only iterated to ensure an accuracy of 1 mm y −1 in flow rates.

Discussion
We quantify here for the first time the magnitude of betweencell groundwater lateral flow globally as a function of grid spacing. We find that the global mean lateral flow at the model grid resolution and its scaling with aggregation is reasonably well described by a simple statistical framework for groundwater flow. We find that at grid spacings smaller than 0.1 • (∼10 km), between-cell groundwater flow is comparable in magnitude to recharge, so clearly needs to be represented if water budgets of individual grid cells are to be realistic. Grid spacings this small are increasingly employed in regional climate change studies that seek to resolve topography, sea breezes, and urban areas [31][32][33]. A transition to generally small between-cell flows takes place between 0.1 • (42% of land area with lateral flow a significant component of the groundwater budget by our conservative combined criterion) and 1 • (1.5% of land area). This justifies to some extent the neglect of lateral groundwater flow in current global climate models, which are often run at resolutions coarser than 1 • . As model resolution improves and as regional and mesoscale models become more widely employed to represent specific weather and climate events, the importance of accounting for lateral groundwater flow increases. The significance of lateral groundwater flow could be estimated for the specific resolution and region of interest using numerical simulations such as ours, or more crudely by using the statistical model developed here, in order to motivate decisions about what kind of groundwater representation to include for a particular climate simulation.
Even where lateral water fluxes per unit area are small when summed over a global model grid cell, an appreciable portion of the grid cell may include groundwater-fed wetlands, riparian zones, or oases, which should ideally be reflected in the modeled cell-aggregated soil and vegetation status and land-atmosphere fluxes through some sort of subgrid parameterization [34], such as those proposed to account for variability in soil moisture [35][36][37] and water table [38]. It has been estimated that shallow groundwater may provide a supplemental water source for aquatic systems and vegetation over 22%-32% of the global land area [2], with as much as 50% of global land area simulated to have mean water table depth shallower than 5 m [39]. Expressions for groundwater discharge area fractions as a function of water table depth, topographic relief, and other parameters may be derived from high-resolution modeling such as ours, combined with observational datasets on wetlands and other groundwaterdependent ecosystems.
While the numerical model we used to simulate groundwater flow is state-of-the-art at the global level, it is subject to limitations which make the quantitative estimates of flow patterns given here tentative. (1) The wide variations known to exist in aquifer permeability and conductivity are mostly not represented, in favor of using 'typical' values that give fairly good results when calibrated against a large dataset of water table depths [2]. [40] have mapped near-surface permeability based on typical values determined for lithologies combined with geologic maps; such datasets of geologic variability could be used to produce more refined estimates of groundwater flow. (2) The 30 resolution does not resolve small discharge features, such as most springs and stream channels. Hence, we may overestimate the size of smaller discharge zones, which are represented as an individual 30 pixel when they should in fact be some fraction of one pixel. If the typical length scale of discharge features (s in equation (6)) is thus overestimated, our statistical model implies that mean |L| at a given grid scale may also be overestimated. (3) On the other hand, our model may underestimate discharge zone extent by failing to explicitly consider plant uptake and capillary rise of shallow groundwater [41,42]. (4) Our model ignores time variation in recharge due to climate variability and change. Variation between wet and dry periods in discharge rate D, discharge area fraction f , and lateral flow pathways can be studied using the 'recession curve' of streamflow in different-order streams within a region [43,44]. (5) We also do not account for uncertainty in the mean recharge estimates of [23] (which could be assessed by repeating this analysis with alternative recharge data sets as drivers); anthropogenic disturbance (such as discharge due to groundwater pumping and recharge due to irrigation return flows) [45,11,46]; or the possibility of nested multiscale flow paths from recharge to discharge areas [47,48].
Despite these limitations, the results presented form a starting point for parameterizing coarser resolution models and for comparing with more capable, finer resolution groundwater models in the future. Comparisons with observational evidence presented in [2] suggest that the modeled areas and magnitudes of groundwater convergence (positive L) are broadly realistic in that shallow water tables and ecologically important wetlands are reproduced.
Our idealized statistical model also ignores these important features of real groundwater flows, and additionally neglects variations in recharge over regional to global domains and nonrandom orientation of sink areas (e.g. along stream channels). Nevertheless, it may provide an easily understood starting point for quantitatively considering the role of groundwater flow at different model resolutions and serve as a benchmark for comparing different numerical results in specific landscapes.
Groundwater hydrologists have developed a nondimensional 'water table ratio' based on mean recharge and aquifer geometry and permeability. A high water table ratio denotes 'topography-dominated' conditions where recharge is high relative to permeability and flow is more local, while a low water table ratio denotes 'recharge-dominated' conditions under which groundwater flow between subbasins is more significant relative to recharge [48,17]. Areas in the United States where large interbasin groundwater transfers were found [21] have been shown to be those with recharge-dominated conditions [49]. It will be interesting to compare the scaling of lateral flow for areas with different values of this water table ratio across climate zones and topographic regimes.

Conclusions
We present a global estimate of the magnitude of between-cell groundwater flow as a function of grid cell size by aggregating results from a model of equilibrium groundwater flow run and validated globally. While our model results are subject to some limitations, they suggest that between-cell groundwater flow needs to be represented in long-term integrations of higherresolution regional models to ensure realistic water budgets, but will have small effects on water exchanges in current global climate models. For coarse grid resolutions, accounting for within-cell groundwater flow through parameterization of subgrid moisture heterogeneity would likely be important.