Simulating the sensitivity of evapotranspiration and streamflow to large-scale groundwater depletion

See allHide authors and affiliations

Science Advances  19 Jun 2019:
Vol. 5, no. 6, eaav4574
DOI: 10.1126/sciadv.aav4574


Groundwater pumping has caused marked aquifer storage declines over the past century. In addition to threatening the viability of groundwater-dependent economic activities, storage losses reshape the hydrologic landscape, shifting groundwater surface water exchanges and surface water availability. A more comprehensive understanding of modern groundwater-depleted systems is needed as we strive for improved simulations and more efficient water resources management. Here, we begin to address this gap by evaluating the impact of 100 years of groundwater declines across the continental United States on simulated watershed behavior. Subsurface storage losses reverberate throughout hydrologic systems, decreasing streamflow and evapotranspiration. Evapotranspiration declines are focused in water-limited periods and shallow groundwater regions. Streamflow losses are widespread and intensify along drainage networks, often occurring far from the point of groundwater abstraction. Our integrated approach illustrates the sensitivity of land surface simulations to groundwater storage levels and a path toward evaluating these connections in large-scale models.


Human development has markedly altered hydrologic systems across the globe. Storing and diverting large volumes of water to support human use and crop production redistributes water spatially, changes streamflow timing, and decreases overall water availability. Recent global studies highlight the extent of human intervention in hydrologic systems and demonstrate large water demand with the expansion of irrigated agriculture (1), as well as hydrologic disconnection caused by river impoundments (2) and cascading water scarcity with human activities (3). Despite local and global efforts to quantify human diversions, many interdependencies between anthropogenic water use and natural hydrologic processes remain poorly quantified. Uncertainty in these connections limits our ability to sustainably manage water resources.

Groundwater is a critical supply for human systems. It supports more than 40% of irrigation globally (4) and is the sole water source for crops in many arid regions that could not otherwise support agriculture. When used in conjunction with surface water, it can stabilize total water supply. Unfortunately, the reliability of groundwater, and its relative abundance in otherwise water-limited locations, has also lead to its overuse. Sustained groundwater pumping has resulted in sustained groundwater storage losses globally (5, 6). In the United States alone, roughly 800 km3 of water was depleted from groundwater storage over the 20th century (7). Unlike surface water reservoirs and streamflow, which are routinely subject to large interannual fluctuations, groundwater storage losses have generally been increasing over time. While multiple studies have evaluated the ability of remaining groundwater supplies to sustain future human demands, outside of several heavily studies regions, much less is known about the impacts of widespread storage losses on watershed function, and global models generally do not include these losses in their simulations.

Groundwater pumping is not an isolated activity; in managed agricultural systems, pumping generally supports irrigation. Connections between irrigation, recharge, soil moisture, atmospheric water content, and downwind precipitation have been demonstrated [e.g., (810)]. However, the impact of storage changes has not been previously isolated from other water management operations (e.g., irrigation and surface water diversions) at the continental scale. Groundwater storage changes are different from many other management impacts because the effects of storage losses are spatially diffuse and temporally persistent; pumping impacts can be observed outside of areas where groundwater irrigation is directly applied and will persist even if irrigation practices change and pumping is curtailed. The long-term storage losses caused by a century of groundwater development can be viewed as a large-scale reshaping of the integrated hydrologic landscape. This will influence watershed response to both natural and human perturbations moving forward. This study seeks to isolate the impact of decreased groundwater storage on the hydrologic landscape and start to unravel the hydrologic differences between modern depleted groundwater systems and natural watersheds.

Groundwater surface water exchanges play an important role in the dynamics of natural hydrologic systems, and there is an increasing push to include groundwater processes in global-scale land surface and climate models. In an idealized sense, groundwater can be conceptualized as a subdued replica of topography, with recharge occurring across the landscape and lateral flow in the subsurface creating convergence zones at low elevations [e.g., (11)] (Fig. 1A). Where the water table is deep (often at high elevations), surface water generally recharges down to the water table; however, when groundwater is within a critical depth range (less than about 10 m from the land surface), connections between water table depth and soil moisture can help support evapotranspiration (ET) (1214). In addition, water can move laterally through the subsurface, from the point of recharge eventually converging to surface water bodies where it can discharge to streams as baseflow [e.g., (15, 16)]. The relatively slow processes of infiltration and lateral redistribution via groundwater flow result in groundwater levels, and groundwater surface water exchanges, which generally respond more slowly and less markedly to surface water changes than other system components (17, 18). As a result of this “long memory” characteristic, groundwater contributions can help stabilize ecosystems by maintaining streamflow and supporting plant productivity during dry periods.

Fig. 1 Conceptual model of groundwater surface water interactions along a hillslope for two groundwater configurations before and after groundwater depletion.

(A and B) The solid black line shows the land surface for an idealized transect. Blue shading indicates saturated groundwater in the subsurface and ponded water (i.e., streamflow) at the surface. The dashed line identifies the depth above which ET is correlated to water table depth. This is illustrated here as a constant depth below the land surface but in a real system could vary based on land cover type and subsurface properties. When the water table (top of the blue shading) is above the dashed line, changes in the water table depth affect soil moisture and water availability for ET (labeled as Groundwater connection in the figure). Where there is ponded water, and the groundwater is directly connected to the surface, we show exchanges driven by the gradient between the stream and the adjacent groundwater levels (red triangle). In the depleted groundwater scenario, the water table is drawn down (dark brown shading) and groundwater surface water connections are affected accordingly.

In developed systems with a history of groundwater pumping, deeper water tables can alter lateral flow, induce recharge, and decrease discharge from the groundwater system as the water table moves toward a new dynamic equilibrium (19, 20). In addition, lower water tables decrease the area of the landscape where the water table is sufficiently shallow to support ET (Fig. 1B). Sustained water level declines can reduce soil moisture and increase the planetary boundary layer height, both of which affect spatial and temporal patterns in ET (9, 21). Regional studies have also connected groundwater overexploitation to changes in circulation in the lowest levels of the atmosphere [e.g., (22, 23)], streamflow declines [e.g., (24)], increased irrigation demand (21), and future drought risk (25) and have found the impacts of pumping and irrigation on latent heat flux and recharge to be similar to 2.5°C warming (26).

We study the impact of long-term groundwater storage declines on simulated surface water behavior across the United States, using an integrated hydrologic model to isolate groundwater depletions from other anthropogenic activities. Developing high-resolution, large-scale integrated hydrologic models is an identified challenge within the hydrologic community (14, 27). Incorporating groundwater processes into continental scale simulations is particularly challenging because (i) subsurface properties are spatially variable and difficult to observe, (ii) groundwater use is often not monitored or reported directly, (iii) long travel times and uncertainty in groundwater response times make it difficult to predict response, and (iv) groundwater surface water interactions can be highly nonlinear and computationally demanding to simulate. As a result, groundwater is often excluded or heavily parameterized in large-scale studies and the role of storage depletions within large-scale simulations has not been previously quantified. Here, we quantify the impact of groundwater development on simulated groundwater surface water interactions across the continental United States (CONUS).


We apply an integrated hydrologic model that simulates hydrologic processes from the groundwater to the top of the canopy at high spatial resolution to study the widespread impacts of groundwater pumping on streamflow and ET. This physically based approach is computationally demanding but also uniquely suited to evaluate the impacts of groundwater development on the hydrologic cycle because it does not ignore or parameterize groundwater surface water exchanges. We evaluate the impacts of groundwater development using a high-resolution (1 km2) model of the majority of the CONUS (~6 million km2) (28) using two scenarios: a predevelopment case without any human activities and a groundwater development case that includes storage losses from pumping. Other than storage loses, no additional anthropogenic activities (e.g., irrigation or surface water diversions) are added to the depletion scenario. The purpose of this work is not to predict the total impact of anthropogenic activities on hydrologic systems but to isolate the impact of groundwater storage losses from other land use trends to better quantify the incremental changes caused by this development. Complete details of the numerical experiments are provided in Materials and Methods.

The storage difference between the predevelopment and developed scenarios at the end of this initialization period (i.e., the start of the second year of simulation) is shown in Fig. 2A. The total storage loss shown here is roughly 800 km3, which is slightly less than the storage losses in the major aquifers in the domain estimated by Konikow (7) plus the storage losses in the unincorporated areas calculated from depletion rates (6). The largest declines are focused in the High Plains Aquifer, but storage losses are also shown across the western United States, in the Mississippi Embayment and parts of the Ohio River basin. The storage losses shown here are not a prediction of water level changes; rather, they are a reflection of current estimates of storage losses and pumping rates. We use these estimates to generate a realistic starting point to evaluate the potential impacts of large-scale subsurface losses on hydrologic system behavior across a range of hydrologic settings and climates in the United States. Thus, the storage changes shown here should be interpreted as the single perturbation applied to an otherwise natural hydrologic system (refer to Materials and Methods for details on the simulation parameters). We intentionally isolate storage losses from the other human development activities that drive groundwater demand, such as urbanization and irrigation, to uniquely quantify the impact of groundwater development on the hydrologic landscape. In the analysis that follows, we evaluate the sensitivity of simulated ET and streamflow to these simulated storage losses.

Fig. 2 Groundwater storage losses and resulting decreases in streamflow and ET between the predevelopment and groundwater-depleted simulations for the simulated conditions.

(A) Storage losses are the difference in the total subsurface storage between the groundwater-depleted and predevelopment simulations at the start of the transient simulation. As described in Materials and Methods, pumping-induced storage losses were calculated on the basis of previous studies of aquifer depletions and groundwater pumping rates. This storage change map is not intended to be predictive. Rather, this is reflective of the calculated pumping storage losses that were applied to the model. This should be interpreted as the perturbation, which is subsequently used to evaluate sensitivity to groundwater storage changes. Both storage losses and ET are plotted on the 1-km2 grid resolution of the model and are reported in length units to reflect an equivalent ponded water depth (A to C). However, groundwater extraction data were not available at this high resolution; therefore, the smaller-scale variability reflected here is caused by heterogeneity in subsurface properties, topography, and climate. Streamflow changes are mapped to river reaches and calculated as the change in flow at the downstream end of each reach relative to the baseline average flow in cubic meters per second (CMS).

Groundwater declines decrease surface water availability

As would be expected, results show decreased streamflow and ET (Fig. 2, B and C) in response to groundwater abstraction (Fig. 2A). Annual volumetric streamflow declines of ~10 to 50% are routinely simulated across the western portion of the domain. In the High Plains, where storage losses are greatest, there are also many locations with streamflow declines greater than 50% and a number of small tributaries that dry up completely.

As noted above, the single perturbation approach applied here is not designed to recreate historical time series; however, we can still compare qualitatively with historical observations of groundwater pumping impacts. Similar streamflow declines and loss of small tributaries have been observed in multiple historical reconstructions of streamflow across the High Plains (2932). For example, one study found declining streamflow in all of the stations evaluated in Nebraska, 85% of stations in Kansas, 50% of stations in Oklahoma, and 33% of stations in Colorado over the period of record (17). Similarly, an analysis of lake inflows found a 99% decline in inflow to the Optima Lake in the Oklahoma panhandle since the 1960s and 65 to 92% declines in the four westernmost reservoirs in Kansas since the 1950s (31). Another study of the Republican River Basin in Southwestern Nebraska found a 75% decrease in the mean annual flow of the Republican River near its entry to Kansas (30). This previous work has correlated declining streamflow observations with increasing trends in groundwater-supported irrigation. Our study shows changes consistent with previous studies in the High Plains and expands beyond this to other regions where detailed analysis of groundwater-pumping impacts has not been previously completed.

Declines in annual ET are spatially distributed with storage losses and are greatest in the more arid western portion of the domain. While the baseline latent heat flux and transpiration partitioning have been validated with observations in previous work (33), the connection between large-scale groundwater depletion and ET has not been studied in the field because of sparse ET observations. However, our results are consistent with previous studies that have demonstrated connections between water table depth and ET. For example, Szilagyi et al. (13) evaluated relationships between shallow groundwater depths and ET in the North Platte and found a roughly 15% reduction in ET as water table depth increased throughout the critical depth range.

It should also be noted that ET losses in this simulated framework occur because we have intentionally isolated groundwater storage changes from all other activities (i.e., irrigation is not being simulated). In reality, the storage losses are generally caused by pumping to support irrigated agriculture. In these systems, irrigation is a dominant control on ET, working to enhance natural plant water availability (8). Therefore, a decrease in summertime peak ET caused by groundwater pumping in isolation, as shown here, could also be viewed as an increased irrigation demand, assuming that water is applied to crops to meet whatever moisture deficit is not supplied naturally (21).

Sensitivity varies with climate and groundwater level

The impacts of groundwater pumping on surface water fluxes (i.e., streamflow and ET) are not necessarily correlated to the magnitude of the storage loss. Sensitivity of surface water fluxes is spatially variable. Annual decreases in total streamflow, local streamflow accretions, and ET are routinely greater than 10% of the long-term storage loss, but the spatial patterns differ depending on which variable you are considering (Fig. 3). While the conceptual model for groundwater surface water interactions is simple, the physical controls of groundwater configuration and water and energy partitioning at the land surface are complex and vary with climate, topography, land cover, soil, and geology (15, 21, 33, 34). The advantage of a large-scale high-resolution model is that we can evaluate the impact of storage changes across a range of climates and hydrologic settings.

Fig. 3 Total surface water declines over the 1-year simulation period relative to long-term storage losses.

Maps of the change in (A) total surface outflow, (B) surface flow accretions, and (C) annual ET between the depleted and predevelopment scenarios, as a fraction of the initial storage loss in the depleted scenario (Fig. 1A). Results are aggregated into roughly 25,000 subbasins (white outlines) that each encompass a single stream reach. Outflow is the total streamflow exiting each subbasin and is therefore an aggregated measure of all of the inflow upstream of the basin. Accretions are the difference between outflow and inflow across a reach (i.e., the stream gains or losses within a subbasin).

Impacts to total streamflow (i.e., outflow at the end of each reach; Fig. 3A) necessarily concentrate along the channel network. Therefore, the stream response is not strongly correlated to, and can far exceed, local storage changes (as indicated by the fractions greater than one, where overland flow losses exceed local pumping). This is caused by upstream streamflow declines propagating through the surface water network as well as lateral groundwater flow. Groundwater surface water exchanges along surface water bodies are a reflection of both local and regional groundwater gradients that control lateral groundwater flow (Fig. 1). Lateral redistribution of water in the subsurface means that water bodies can be affected by groundwater storage changes that are far away, as long as the system is hydrologically connected.

The role of groundwater surface water exchanges in streamflow changes is further illustrated by the fractional changes in streamflow accretions (Fig. 3B). Local streamflow accretions are defined as the net gains and losses along a given stream reach (i.e., the difference in flow from the beginning of a reach to the end). Accretions can result from local precipitation entering the stream or exchanges with groundwater as either baseflow gains or recharge losses. The largest impact to streamflow accretions occurs in the eastern portion of the domain, where shallow groundwater levels in the predevelopment simulation result in greater connection between groundwater and surface water bodies and more topographically driven groundwater flow (34). The relative impact of storage losses on streamflow accretion is smaller in the western portion of the domain, but changes in annual streamflow accretions still routinely exceed 10% of the local long-term storage loss. This illustrates how a single volumetric loss to the system can be perpetuated through recurring annual changes at the surface as groundwater surface water interactions are shifted.

In contrast to streamflow changes, ET losses (Fig. 3C) are spatially diffuse because ET varies with the local water table depth (Fig. 1). The sensitivity of ET to storage changes is spatially variable but has the opposite pattern to streamflow accretions; sensitivity to storage changes is generally greater in the more arid western portion of the domain, which is already water limited in all or part of the year. In addition, when the simulated natural water table depths are less than roughly 5 m, there is a strong positive correlation between the total storage change in a subbasin and the overall ET decrease. This is consistent with previous studies that illustrate a dominant groundwater control on ET sensitivity in this depth range [e.g., (8, 35)].

Groundwater recharge increases and discharge decreases

In addition to spatial variability, groundwater surface water exchanges also have strong seasonal oscillations. In general, groundwater storage increases during wet periods as excess moisture is recharged to the subsurface and decreases during dry periods as groundwater is discharged to streams as baseflow and depleted by ET. Monthly fluctuations in groundwater storage across the domain can be as much as 150 km3 (Fig. 4D); while this is large relative to the long-term storage loss of 800 km3, it should be noted that winter recharge is balanced by summer discharge and the annual trend is not this high. The exact timing of seasonal oscillations varies spatially but aggregated across the domain; groundwater recharge peaks during the winter and discharge peaks during the summer. This is consistent with maximum ET and streamflow, which occur in the summer and late spring, respectively.

Fig. 4 Seasonal decreases in net groundwater contributions to the land surface.

Shaded areas in the maps (A to C) show locations where the groundwater development scenario results in a smaller net contribution of groundwater to the surface water system relative to the predevelopment scenario. Contributions are calculated and plotted on a grid cell basis using the groundwater storage changes over each season. Decreased contributions can occur because of either increased groundwater recharge (purple), decreased discharge (green), or a swap of discharge to recharge (blue). (D) Total monthly ET, streamflow, and changes in storage summed across the domain for the predevelopment simulation. Arrows indicate the overall direction that groundwater depletions shift each component. Note that groundwater pumping is not included in (D) because there is no pumping in the predevelopment simulation, but there is an additional flux of 1.6 km3 out of the domain per month out to reflect ongoing pumping activities included in the developed scenario.

Connections between groundwater pumping and stream changes are well established in the literature through the concept of “pumping capture.” Before human development, groundwater systems were generally in a state of approximate dynamic equilibrium, where groundwater recharge was balanced by discharge over long time periods (20). Groundwater pumping represents an increased discharge from the system. Continuity dictates that this perturbation must be balanced by some combination of increased groundwater recharge, decreased discharge (also referred to as captured discharge), or changes in groundwater storage (20). When pumping is first initiated, the losses are realized as changes in storage; however, over time, as pumping continues, the system will move toward a new dynamic equilibrium in which extractions are largely compensated for by induced recharge or decreased discharge (19, 20).

Similar to real-world watersheds, our model simulates a semiclosed system driven by meteorology. To maintain continuity within the model, groundwater storage changes must be balanced by groundwater pumping, surface outflows, and ET. This work is unique from many previous pumping analyses because we directly incorporate physically based groundwater connections to ET in addition to streamflow. Therefore, a net “discharge” from the groundwater system can be in the form of not only baseflow to a stream but also ET, which is supported by vertical water movement in the subsurface. In the simulations presented here, the perturbed system has not been run to a new dynamic equilibrium (and it can be argued that developed groundwater systems are not in a state of equilibrium). Rather, we present the changes in recharge and discharge induced by the storage perturbation.

The ET and streamflow declines (Fig. 3) are caused by a combination of increased groundwater recharge, decreased discharge, and locations that were previously net dischargers swapping to net recharge. Increased recharge is most dominant in the winter (October to January; Fig. 4A) when the system is most energy limited. In the western portion of the domain, especially at high elevations, patterns of increased recharge also extend into the spring (February to June; Fig. 4B). Decreased discharge is prevalent in the water-limited summer months across the domain (June to September; Fig. 4C) and in the spring for the lower elevation and more humid eastern portions of the domain. Areas where groundwater discharge swaps to recharge are concentrated in the High Plains and are especially prevalent in the spring. This pattern is consistent with the previous discussion of streamflow declines and loss of small tributaries in this portion of the domain. For example, Kustu et al. (17) demonstrated decreasing trends in annual and dry-season streamflow and an increase in the number of low flow days across the Northern High Plains from 1950 to 1980, the time period that corresponds to intensive irrigation development.

Seasonal variability in the surface system is dampened

The increase in groundwater recharge during the winter and decrease in discharge in the summer caused by storage losses decrease the magnitude of the seasonal oscillation of both ET and streamflow (Fig. 5). Surface sensitivity to groundwater storage losses is greatest during water-limited periods when groundwater supports streamflow and ET (Fig. 4D). Seasonal dampening is primarily caused by decreases in the peak discharge and ET that occur in the spring and summer, respectively. Streamflow and ET are also decreased in the winter as groundwater recharge increases; however, ET is less sensitive to groundwater changes in the winter because the system is generally more energy limited.

Fig. 5 Changes in the seasonal amplitude of streamflow and ET.

(A and B) Maps of the difference in the seasonal amplitude (monthly maximum − monthly minimum) in the groundwater-developed simulations compared to predevelopment simulations. Results are aggregated into roughly 25,000 subbasins that each encompass a single stream reach. Streamflow is the total flow at the outlet, and ET is the total across the subbasin area.

The seasonal amplitude dampening is routinely greater than 10% for both streamflow (Fig. 5A) and ET (Fig. 5B) especially in the more arid western portion of the domain. Seasonal streamflow amplitude declines greater than 25% are consistent with areas where streamflow is markedly reduced (Fig. 2B). Perhaps more significantly though, amplitude declines of 5 to 25% persist downstream along major rivers such as the Missouri, the Ohio, and the Arkansas as changes propagate through the stream network. Similarly, overall losses and seasonal amplitude changes in ET are spatially correlated with each other. In the more humid, eastern portion of the domain, the simulated ET losses (Fig. 2) are small relative to the seasonal variability, and therefore, no amplitude changes are noted despite the fact that storage changes do occur.


Large-scale groundwater depletions are well documented in the United States and across the globe. In addition to threatening the sustainability of groundwater-dependent systems, sustained declines also change surface water availability and watershed function. Groundwater plays an important role as a long memory stabilizer of both human and natural systems. Storage losses reconfigure groundwater resources and can systematically shift groundwater surface water exchanges as watersheds transition to a new dynamic equilibrium with lower storage levels. The simulated sensitivities to groundwater depletion shown here are consistent with local impacts demonstrated by previous work from the groundwater community. However, the large-scale impacts of groundwater development on watershed function and simulated behaviors have not been well studied outside focused watershed analyses. Lateral groundwater flow is often not well represented in large-scale modeling approaches, much less the sensitivity to groundwater development and storage losses. These connections are increasingly important as we try to quantify large-scale trends in hydrology and disentangle impacts of climate change from human water use. For example, one recent study showed that the 100th meridian (the longitudinal band across the Great Plains of the United States where precipitation roughly equals potential ET) has been shifting eastward over the last 25 years as a result of decreasing precipitation and increasing potential ET (36). Here, we show that large-scale storage losses can also contribute to increased aridity and regional decreases in surface water availability.

We perform a numerical experiment to isolate the impact of groundwater storage changes from other human developments. Our analysis quantifies how this single human activity propagates through otherwise undeveloped surface water systems across the CONUS. Historically, groundwater development has occurred in conjunction with other agricultural activities (e.g., land cover change and irrigation), and therefore, our analysis is not predictive of long-term trends in irrigated agriculture across the United States. Rather, we isolate the large-scale anthropogenic reconfiguration of the subsurface from other activities to better understand how today’s groundwater-depleted systems might behave differently than the undeveloped systems that were common a century ago. Evaluating storage changes separately from other anthropogenic activities is important because groundwater impacts extend spatially beyond the point of extraction and persist temporally even if pumping is stopped.

Our simulations show that groundwater storage losses decrease overall surface water availability. Streamflow declines and loss of small headwater streams occur near pumping but also persist downstream. Lateral groundwater flow translates storage losses to stream impacts far from the point of pumping. Evapotranspiration also declines with groundwater losses as a result of deeper water tables, which decrease the areas of strong coupling between groundwater, soil moisture, and ET. Because irrigation is not simulated here, decreased ET can also be conceptualized as a decrease in the naturally supported ET and potentially an increase in irrigation demand. These results illustrate the importance of considering the impacts of groundwater storage levels on integrated system behavior. We show that groundwater depletions alter the exchanges and partitioning that occur throughout the system: increasing recharge, decreasing discharge, and dampening the seasonal variability of both streamflow and ET.

In combination, these results illustrate the role of groundwater development in reshaping the hydrologic landscape. Given the ubiquity of groundwater pumping across the globe, and a growing focus on integrated water management, it is critical that hydrologists work to understand the behavior of developed systems with the same physical rigor applied to natural ones. As these lateral flows are not represented in earth system simulation models, it becomes important for future work to explore how more simplified approaches (e.g., land surface models or water management models) compare to integrated hydrologic simulations. Here, we show that watershed response to storage changes varies with season, climate, and hydrologic setting. We demonstrate one approach to evaluate the propagation of groundwater storage changes through hydrologic systems and highlight the sensitivity of simulated large-scale land surface behavior to subsurface storage levels.


Integrated hydrologic model

Our simulations used the fully integrated physical hydrology model ParFlow-CLM (3739). ParFlow is a finite difference model that solves three-dimensional variably saturated flow in the subsurface using Richards’ equation. Overland flow is simulated using Manning’s equation and the kinematic wave formulation using a finite-volume approach. Groundwater flow is fully integrated with surface water flow using a free-surface overland flow boundary condition. Land surface processes such as plant water use, bare soil evaporation, snow accumulation and melt, radiation and land-energy budget, and soil temperature are simulated using a variation of the Common Land Model (CLM) coupled to ParFlow (12). The combined ParFlow-CLM model solves the full water-energy balance at the land surface and captures physically based hydrologic interactions from the bedrock through the top of the canopy. ParFlow is an established open-source platform that has been well verified and applied in many watersheds across the globe.

While the integrated approach is computationally expensive, it has several advantages. First, ParFlow-CLM does not require any a priori specification of inundated areas (i.e., rivers, lakes, and wetlands). This is because the groundwater and surface water systems are solved within a fully integrated framework; natural surface water bodies and unsaturated zones can expand and contract only on the basis of the fluxes through the system and physical properties such as topography and hydraulic conductivity. In addition, because all water fluxes and stores are derived from the solution of coupled partial differential equations, there is no need to parameterize the exchanges between the surface and subsurface (i.e., with a stream package type formulation). This is also true for the partitioning of precipitation fluxes at the land surface into overland flow and recharge. Recharge fractions are not specified; rather, runoff and recharge are a product of soil properties, moisture content, and topography. Combined, these factors allow ParFlow-CLM to simulate changes in storage, recharge, ET, and discharge that evolve dynamically in response to groundwater storage changes. These capabilities are consistent with other integrated hydrologic models that have generally been applied at the watershed scale but unique compared to typical large-scale modeling approaches such as global water balance models and earth system models, which simplify these interactions and do include lateral groundwater flow.

Numerical experiments

Two numerical experiments were conducted to quantify the impacts of groundwater pumping on streamflow and ET: a baseline simulation with no pumping and a groundwater development scenario with pumping. Both cases were based on the CONUS-ParFlow model (28), which encompasses approximately 6.3 million km2 of North America, focusing on the United States. These simulations have a lateral resolution of 1 km and five subsurface layers ranging in thickness from 0.1 to 100 m for a total thickness of 102 m using a terrain-following grid (38). Subsurface units are laterally homogeneous within each layer. Model inputs were developed using publicly available national data products of soils and geology.

Because groundwater flow dynamics have a long memory, both numerical experiments have a similar progression: A simulation with constant recharge spanning multiple centuries (28) was used to initialize a series of transient simulations driven by hourly historical reconstructed meteorology for an average water year (33). The baseline case has no anthropogenic impacts and thus represents predevelopment conditions; the groundwater development scenario included historical depletions after the initialization step and a modern reconstruction of groundwater pumping (without irrigation) applied during all the transient simulations. Thus, the two simulations differ only by the impacts of groundwater pumping over the last century.

A detailed validation of the predevelopment simulation is provided in the supplementary information of Maxwell and Condon (33). We show that despite remaining uncertainties and data limitations, there is good agreement between simulated outputs and observations. We also show that the physical interactions between groundwater and surface water are appropriately captured and that this approach can be used to better quantify groundwater contribution to land energy fluxes over large scales (33).

For the groundwater pumping simulations, we combined the best available estimates of aggregated aquifer storage losses (7) and spatially distributed pumping rates (6). Roughly 1000 km3 was lost from storage across 40 different aquifers within the United States from 1900 to 2008 (7). Some of the most marked losses occurred in the High Plains and Central Valley as a result of intensively irrigated agriculture (40); however, widespread pumping also occurred outside these areas. Unfortunately, spatially distributed estimates of groundwater depletions are not available over much of the United States; therefore, we combine lumped aquifer estimates of long-term storage losses with spatially distributed groundwater pumping estimates to derive a spatially distributed storage loss estimate.

Konikow (7) calculated lumped storage changes for 40 separate aquifers using a variety of methods including observed changes in groundwater levels and heads combined with estimates of aquifer properties, storage loss estimates from the GRACE satellite data, flow models calibrated to observed heads, water budget calculations based on pumping data and other fluxes, pumping data combined with consumptive use fractions, and volume estimates based on subsidence data. While these storage changes (7) represent a thorough long-term accounting of a total volume of water removed from groundwater, they are not directly usable within the model because the estimated storage losses are aggregated by a major aquifer unit and do not account for pumping outside these aquifers.

Therefore, we started from a spatially distributed groundwater pumping dataset (6) and adjusted the pumping rates to match this total lumped storage loss estimates (7). The groundwater pumping dataset developed by Wada et al. (6) estimated monthly groundwater extractions from 1960 to 2010 at a 6-min (~11 km) lateral resolution for the United States, and 0.5° resolution globally. Their extraction estimates are based on data from the International Groundwater Resource Assessment Centre and average extractions aggregated by county in the United States. Wada et al. (6) downscaled these data to the 6-min grid using estimates of total water demand, including water consumption for irrigation, livestock, domestic and industrial uses, and estimates of surface water availability. This dataset provides spatially distributed information lacking from the aquifer storage change estimates but also is limited because it estimates total groundwater extractions, not storage changes. That is, the extraction data reflect total groundwater pumping but do not account for changes in recharge, so they are not an estimate of the net groundwater depletion that results from pumping.

Here, we used the spatial patterns in average extraction rates from 1960 to 2010 (6) to spatially disaggregate the 1900–2008 aquifer storage losses (7). This approach ensures that the total depletion magnitudes over the aquifer areas are equivalent to the storage losses. For every aquifer, we calculated a depletion factor by summing the total extractions, assuming that the long-term average extraction rate is applied for 109 years and dividing by total storage losses. For grid cells outside the 40 aquifer units, the storage loss was calculated by scaling the extraction estimates from each grid cell by the local depletion fraction, where the local depletion fraction is assumed to be the average depletion fraction for all of the aquifers within the major river basin. The result of these steps is a storage change map, where the spatial distribution of depletions matches the extraction rates (6) and the total volume of the depletions matches the 1900–2008 storage losses (7).

For the groundwater development simulation, the groundwater depletions developed by the downscaling approach described here were applied after initialization, and a constant groundwater pumping flux was also applied during the transient simulation. The groundwater extraction flux was calculated using the gridded extraction data from Wada et al. (6), which represents the average extraction rate from 1960 to 2010.

In addition to the uncertainty and data limitations of the two input datasets, we also note two additional assumptions that were imposed with this approach: (i) The spatial distribution of groundwater extraction does not change over the development period, and (ii) constant depletion fractions for aquifers are representative of outlying areas. As noted previously, the intent of the groundwater development simulations is not to predict groundwater levels in the United States. Rather, the depletions were incorporated into the simulation to create a realistic perturbation to the groundwater system that can be implemented in a systematic way within a numerical modeling framework.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: The simulations presented here were completed using resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under contract no. DE-AC02-05CH11231. Funding: This work was supported by the U.S. Department of Energy Office of Science, Offices of Advanced Scientific Computing Research and Biological and Environmental Sciences IDEAS project and the Sustainable Systems Scientific Focus Area under award number DE-AC02-05CH11231. Author contributions: L.E.C. completed the pumping simulations and prepared the manuscript. R.M.M. led the baseline simulations, provided input on pumping simulations, and assisted with manuscript preparation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. ParFlow is freely available through GitHub (, and additional model resources can be found at The model outputs used for this work are published with the following doi: Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article