The Multi‐Scale Dynamics of Groundwater Depletion

Unsustainable depletion of aquifer storage and diversion of groundwater from downstream users continue to be serious problems, even though their adverse effects are widely recognized. Groundwater depletion involves interactions between economically motivated pumping decisions and physical constraints. Here, we investigate these interactions by using optimal control techniques to describe the pumping decisions of economic agents who share an aquifer. Our approach relies on a multi‐scale description of unconfined groundwater flow, applied to a computational experiment that illustrates some important impacts of aquifer development. We show that cooperative groundwater management can provide higher economic benefits, with less storage depletion, than an uncooperative alternative. However, demand for water still drives pumping decisions, even when users cooperate. In many aquifers the benefits of pumping are necessarily accompanied by a reduction in aquifer outflow and by depletion of the groundwater reserve, for both cooperative and uncooperative management. These benefits decrease significantly when pumping is limited by well yield constraints or when the diversion of aquifer outflow is restricted. Our overall analysis emphasizes the importance of viewing aquifer management as part of a larger resource allocation problem that considers the conflicting needs of well operators, downstream users, and future consumers of the resource.

analysis range from the diameter of an individual well (typically less than 1 meter) to overall aquifer dimensions (typically 10-1000 km). Relevant time scales range from days to many decades. Our methods accommodate this wide range of time and space scales while recognizing the practical computational limits of nonlinear optimization algorithms.
In this study, we formulate and solve a set of model-constrained optimal control problems that provide quantitative responses to our research questions. We show, with a computational experiment, how this integrated optimization-modeling approach can be used to assess groundwater depletion for different management strategies and aquifer configurations. This experiment is designed to illustrate generic issues in a straightforward way. However, our general approach is sufficiently flexible to deal with variants that include additional physical processes, different agent objectives, or more geological complexity for particular sites.
There is a large body of literature dealing with applications of modeling, optimal control, and game theory methods to the groundwater depletion problem (for example, Bredehoeft & Young, 1970;Brozovic et al., 2010;Burt, 1966Burt, , 1967Gisser & Sanchez, 1980;Madani & Dinar, 2012a, 2012bMüller et al., 2017;Negri, 1989;Provencher & Burt, 1993;Rubio & Casino, 2001, 2003Sahu, 2018a). Much of this literature poses groundwater depletion as a type of common pool problem (Gardner et al., 1990;Hardin, 1968;Madani & Dinar, 2012a;Ostrom, 1990). Bierkens and Wada (2019) summarize some distinctive features of groundwater common pool problem that are relevant to our analysis. In particular, they point out that common pool concepts are most applicable to extensive aquifers that are readily accessed by many agents (e.g., landowners) who are significantly impacted by the pumping of other agents. That is the situation considered here.
In a common pool setting multiple agents may have different objectives and may or may not cooperate in their efforts to optimize their respective benefits. The multi-objective nature of the common pool decision problem significantly complicates the derivation and interpretation of uncooperative behavior (Başar & Olsder, 1998;Fudenberg & Tirole, 1991). It is helpful to simplify the problem by bounding the range of possible pumping strategies with respect to two factors: cooperation and foresight.
On one extreme is a strategy that relies on complete cooperation and perfect foresight. Such a strategy gives pumping decisions that maximize the total net revenue obtained by all agents over a multi-year operating period. These agents are assumed to have perfect knowledge of the future consequences of their current pumping decisions. The benefits obtained through cooperative operation are distributed among the agents by mutual consent. In effect, the aquifer is managed by a single meta-agent who seeks to optimize aggregate benefit.
On the other extreme is a strategy that relies on many uncooperative agents making independent decisions with no foresight. This strategy, which is sometimes termed "myopic," gives pumping decisions that separately maximize the current benefit for each individual agent, without regard for future consequences and without any consideration of the actions of other agents. This process is repeated at all decision times in the operating period. For simplicity, we use the term "cooperative" to refer to the cooperative-with-foresight strategy and "uncooperative" to refer to the uncooperative-without-foresight strategy. Brozovic et al. (2010), Negri (1989), Provencher and Burt (1993), and Rubio and Casino (2001) argue that a cooperative decision strategy gives greater total net revenue from a groundwater resource than the uncooperative alternative, but Gisser and Sanchez (1980) present an example that shows that the difference between the two options is negligible. All of these studies rely on highly simplified aquifer models that cannot adequately portray the benefits and limitations of cooperation because they do not consider critical physical processes. The more realistic multi-scale aquifer model used here enables us to consider factors such as the impact of local drawdowns at individual wells, gradual declines in regional water tables, capture of recharge, well yield limitations, intermittent pumping, the nonlinear dependence of transmissivity on unconfined water levels, and the effect of sustainability constraints that limit pumping. These model capabilities provide a clearer picture of the effects of management strategy and aquifer characteristics on sustainability metrics as well as revenue. The following section describes our general approach and summarizes the computational experiment we use to investigate our research questions.

Methodology
An optimal control approach for analyzing groundwater depletion requires specification of agent objective functions as well as dynamic constraints that restrict the set of feasible solutions. The unknown decision variables in this problem are annual pumping rates that are determined in real time at multiple wells scattered throughout an unconfined aquifer. Each well is operated by a different agent over a 25 to 50-year time period. The wells are pumped continuously at a constant rate during each six-month growing season, with no pumping during the remaining 6 months of the year. This assumption illustrates the effect of seasonally variable pumping and can easily be modified to deal with other schedules.
The pumping rates and their impacts on aquifer storage and outflow are determined by the optimization procedure. The objectives we use for the cooperative and uncooperative strategies differ with respect to (i) the procedure for allocating benefits to individual agents and (ii) the information available for making a decision. The optimization constraints are derived from our multi-scale groundwater model and from optional conditions imposed on aquifer outflows. Mathematical summaries of the model and the optimization procedure are provided in the Supplementary Information (SI).
We use our optimization procedure to investigate aquifer management strategies for several different physical situations. The annual components of the multi-agent cooperative objective and of the individual agent uncooperative objectives all express net revenue as the difference between the demand for pumped groundwater and the cost of extracting this water. The demand is a concave quadratic function of pumping rate and the cost is a bilinear function that depends on the product of pumping rate and well lift (Negri, 1989).
We evaluate the relative performance of the two managemnt strategies by comparing the discounted cooperative objective value over the entire operating period with the discounted sum of the annual uncooperative objective values, summed over all agents. These are commensurate measures of the net present value benefit extracted from the resource. We supplement this revenue-based performance measure with a number of other important sustainability metrics. These include the cumulative volumes of water pumped, storage depleted, and recharge captured. Taken together, these metrics describe the dynamics of the depletion process and quantify tradeoffs between net revenue and environmental impacts.
The aquifer model deserves some elaboration since it provides the multi-scale capabilities we need to properly describe the depletion process. This model must be able to simulate short-term drawdowns at intermittently pumped small-diameter wells, since these affect well yield and pumping cost through the well lift. It must also accurately describe the long-term regional head distribution, which affects local well lifts and controls processes such as the capture of aquifer outflow. We deal with this wide range of scales by merging local head solutions that apply near pumping wells with a regional head solution that applies away from these wells. The general approach is summarized in Figure 1.
The local head near each well is obtained from a numerical solution of a one-dimensional axially symmetric unconfined groundwater flow equation that describes vertically averaged flow near a fully completed well (Guymon, 1980). This equation relies on the Dupuit-Forchheimer assumption that vertical flow resistance can be ignored (Kirkham, 1967). Haitjema (2006) provides guidelines that indicate this assumption is justified in our application, where the wells are fully completed and far from aquifer boundaries. The well equation is nonlinear when flow is unconfined since the transmissivity depends on the head through the saturated depth.
The near-well groundwater flow equation is solved on a daily time scale for each well inside a circle of specified radius from the well center, as shown in Figure 1b. In our computational example, the numerical discretization along the line from the well center to the circle boundary uses 50 grid points that are spaced logarithmically from 0.5 (the well borehole radius) to about 7 m (the largest grid spacing for the local solution). The concentric circles shown inside the well circle pass through the near-field radial grid points. This discretization approach provides a gradual transition from the small scale of the borehole to the much larger scale of the regional grid (discussed below).
The head boundary condition on the well circle is obtained from a finite element method (FEM) numerical solution of the two-dimensional regional flow equation (Pinder & Gray, 2013), which is discretized over the entire aquifer on a polygonal spatial grid, as shown in Figure 1a. This equation also relies on the Dupuit-Forchheimer assumption and is solved at a daily time scale. The nodes on the triangular FEM grid that are closest to any given well are arranged to lie on the corresponding well circle boundary. The individual near-well solutions apply inside their respective well circles while the regional FEM solution applies everywhere outside the circles.
The transition from the axially symmetric near-well head solution to the asymmetric regional solution is an approximation that needs to satisfy mass balance conditions on each well circle, to a reasonable degree of accuracy. Water mass is conserved across the well circle boundary separating the near-well and regional grids if the head and head gradient are continuous at every point on this boundary. In our implementation the near-well (inner) model head on the circle boundary is set equal to the spatial average of the regional (outer) head on this boundary. The outer head boundary values at particular points generally vary somewhat from this average. The continuity of the inner and outer gradients is not guaranteed but can be checked by evaluating mass fluxes across the boundary. Both gradients depend on the well circle radius, which can be adjusted to obtain acceptable mass balance results. The overall accuracy of the head and head gradient approximations at multiple times and locations can be assessed by comparing the annual volume pumped from each well to the sum of the net annual volume flowing inward across the well circle boundary (computed from either inner or outer head solutions) plus the net annual storage volume lost inside this boundary. In our example, this check gives annual volume discrepancies of a few percent for a well radius of 500 m at each of the nine wells.
The near-well and regional groundwater equations used in our multi-scale aquifer model are both nonlinear. This makes the optimization problem nonlinear as well. Other multi-scale groundwater modeling approaches have been proposed by Charbeneau and Street (1979), who consider only linear steady-state problems, and by Foster et al. (2017) and Upton et al. (2019). The approach described here is especially well-suited for integration with readily available nonlinear optimization software. The multi-scale groundwater model can be incorporated into our nonlinear optimization procedure by inserting all of its space and time-discretized equations as equality constraints. This is inefficient since the model equations are solved throughout the aquifer at a daily time step but we only need the annual pumping costs at the wells and the annual minimum aquifer outflow for optimization. Consequently, we use an iterative multi-well version of the response matrix method (Gorelick, 1983;Greenberg, 2015;Maddock, 1972) to derive annual constraints that relate pumping rates directly to well costs and aquifer outflows. The resulting nonlinear optimization problem is conveniently solved with a sequential quadratic programming algorithm (Nocedal & Wright, 2006). Our version of the response matrix approach increases computational efficiency over traditional methods by using a sparse (block-diagonal dominant) approximation of the response matrix within the nonlinear programming iteration loop.
The decrease in regional head that occurs over long times when an unconfined aquifer is depleted gradually reduces transmissivity at the pumping wells by decreasing saturated depth, which we define as the difference between the water table elevation and the aquifer bottom elevation (see Figure 1c). The aquifer and well bottom elevations coincide in our application, which considers only fully completed wells. As the saturated depth and transmissivity decrease at a pumping well there is a gradual reduction in the maximum pumping rate, or well yield, that can be sustained without dewatering the well (Foster et al., 2015;Hecox, 2002;Upton, 2019). Dewatering can be accurately described only if the near-well portion of the multi-scale model can resolve distances on the order of the well radius. We use this capability to impose dynamic pumping rate constraints that prevent dewatering. Additional pumping rate constraints can also be included in the cooperative case to ensure that the total aquifer outflow always remains above a specified minimum value. Well yield and aquifer outflow constraints become economically important when they force pumping rates to take on values below the unconstrained values that maximize net revenue.

Computational Experiment
In our investigation of groundwater depletion, we focus on a computational example that illustrates a number of important features that are often missed in studies based on simpler models. In particular, we examine the effects of management strategy (cooperative vs. uncooperative), economic incentives, well yield limitations, and aquifer outflow constraints. These effects are illustrated with three case studies that rely on three distinct sets of model inputs.
Our example is based on a 15 km by 10 km rectangular unconfined aquifer extending along an east-west axis, illustrated in the plan view in Figure 2. The aquifer is assumed to be in a semi-arid area where groundwater is the major source of irrigation water during a growing season with little or no rainfall and high potential evapotranspiration. Such conditions occur in parts of the western US, as well as in many other regions, and are primary areas of concern for groundwater depletion (Bierkens & Wada, 2019).
Natural long-term average groundwater recharge is introduced through the eastern portion of the northern boundary in order to simulate a water source that originates from rainfall and snowmelt at higher elevation regions outside the valley aquifer. The central portion of the western boundary is controlled by a specified head condition. The edges of this boundary, the western portion of the northern boundary, and the eastern and southern boundaries are assigned no-flux conditions. Aquifer recharge from direct rainfall and irrigation returns is assumed to be small since the aquifer is in a semi-arid area with a relatively deep water table, following discussions by Gee and Hillel (1988) and Jiménez-Martínez et al. (2009). Values for these fluxes, which generally vary over time and space, could be added in situations where they can be estimated and are believed to be comparable to other problem fluxes.
The horizontal hydraulic conductivity is assumed to be uniform throughout the aquifer but the transmissivity that controls two-dimensional flow varies with saturated depth in both the near-well and regional parts of the aquifer. Each of the nine agents operates one of the pumping wells that are scattered throughout the southern portion of the aquifer. In our example, these agents all start pumping at the same time but they could be phased in over the operating period to simulate more gradual development of the aquifer. The well locations and identification numbers are shown in Figure 2.
The aquifer bottom elevation slopes upwards along the x coordinate from the well field to the western outflow boundary, reducing transmissivity in the shallow section of the aquifer that lies near this boundary. Although the agents have the same demand and cost functions, they encounter different transmissivities and well lifts that can cause pumping rates to differ at their respective wells. Table S1 in Supporting Information S1 summarizes information about aquifer dimensions, hydrogeological properties, and other input values used in the computational experiment.
Under unpumped steady-state conditions all of the recharge from the northern boundary converges to the central specified-head portion of the western boundary, where it flows out of the aquifer (Figure 2a). This groundwater outflow could help maintain environmental flows to downstream surface waters, it could be pumped by downstream groundwater or surface water users, and/or it could help prevent salinity intrusion from nearby coastal waters. We use a specified-head outflow boundary as a generic way to account for situations where the aquifer flux leaving this boundary is reduced by pumping within the aquifer, while still remaining positive. The outflow boundary condition could be modified to deal with other situations, such as a boundary formed by a stream that is initially gaining but eventually becomes losing and hydraulically disconnected from the aquifer. In this case inflow from the losing stream to the aquifer (negative aquifer outflow) approaches a constant value independent of the groundwater head at the boundary (Brunner et al., 2011, Figure 1). When this inflow can be estimated such cases can be handled by using a specified flux condition along the central western boundary. The unpumped steady-state head distribution provides the initial condition for an operating period that begins when aquifer pumping starts. When the aquifer is pumped some or all of the northern boundary recharge is captured by the pumping wells. This reduces the outflow leaving the western boundary. Pumping also reduces groundwater storage, as revealed by the drop in water levels, as shown by comparing Figures 2a and 2b. Pumping decisions depend on a two-way iteraction between local and regional water levels. The local well drawdowns that develop over each year's 6-month pumping season recover significantly, but not completely, over the subsequent 6-month unpumped season. This incomplete recovery gradually leads to regional storage depletion and water level decline. Lowering of regional water levels in turn increases pumping lifts and costs at individual wells. This gradual increase in cost can lead agents to decrease their pumping over time.

Results
In this section, we use the modeling and optimization capabilities described above to examine three case studies. In the first case, groundwater pumping and related depletion are determined only by the demand for water and pumping cost. In the second case, pumping is also constrained by well yield, which is related to the aquifer's hydraulic properties, the regional head distribution and the aquifer depth, as described in Text S4. In the final case, pumping is also constrained by an outflow limit designed to protect downstream water uses. The well yield and outflow constraints limit the economic benefits that users can extract from the groundwater stock. The three cases rely on different aquifer properties, economic parameters, and environmental constraints. These inputs have been selected to illustrate different outcomes that might be observed for an aquifer similar to the one described above.
Case 1. Pumping constrained only by water demand and cost Figure 3 provides a convenient overview of the economically constrained case by showing how the daily aquifer fluxes (all measured in m 3 day −1 ) vary over a 50-year operating period. The effect of intermittent pumping is most evident in the seasonal oscillations in the total pumping rate (orange), which is summed over all agents, and in the aquifer storage change (blue). Note that the annual average pumping rate (measured in m 3 day −1 ) is the constant daily flux pumped over the 182-day irrigation season multiplied by the fraction of the year assigned to that season (182/365). The constant daily recharge rate is indicated by the horizontal red line. Daily groundwater storage change (+for storage increases and-for decreases) is computed as the residual in the daily aquifer-scale water balance, with the values confirmed by a mass balance check that considers temporal changes in the regional hydraulic head distribution. Aquifer storage is depleted during the 6-month pumping season and partially restored during the unpumped season. The daily aquifer outflow (maroon, negative for flow out) slowly decreases in magnitude as more water is captured by pumping.
The two columns of plots in Figure 3 compare low-demand and high-demand responses for cooperative and uncooperative strategies. The upper four plots apply to a shallow aquifer (h b = −50 m) while the lower two apply to a deeper (higher transmissivity) aquifer (h b = −100 m). Table 1 summarizes input values and revenue results for the shallow aquifer alternatives.
All of the aquifer pumping rates shown in Figure 3 drop somewhat over time, reflecting the effect of gradually decreasing well heads and associated increases in pumping lifts and costs. When optimum pumping rates fall, the agents must either obtain irrigation water from other sources or modify their cropping strategies (e.g., by changing crops, reducing cultivated area, and/or reducing applied water). The effect of water demand on pumping is substantial (compare the upper and middle rows of Figure 3), reflecting the important influence of economic factors on aquifer depletion. A 50-meter increase in aquifer depth has relatively little impact on the high demand pumping trajectories (compare the upper and lower rows). The uncooperative strategy consistently gives higher pumping rates that not only generate higher total revenue but also incur higher pumping cost, with the overall effect being lower net present value benefit and higher depletion. The net present value benefits are compared in Table 1. Figure 4 shows cumulative totals of the annual pumping, storage depletion, and recharge capture volumes up to the indicated time for high and low demand with a shallow aquifer. The cooperative strategy is shown in solid lines, and the uncooperative strategy is shown in dashed lines. The cumulative capture volume is the cumulative recharge volume minus the cumulative outflow across the western aquifer boundary. Although storage depletion (blue) provides most of the pumped water in all the cases shown here, the contribution of captured recharge (red) increases over time.
These cumulative plots confirm that pumping reduces outflow to downstream users while also reducing the groundwater stock, even in low demand situations where the annual volume pumped is less than annual recharge. This occurs with cooperative as well as uncooperative management. In all the cases shown, recharge is gradually captured by the pumping wells as they generate south-pointing gradients from the recharge boundary to the pumping field.
Pumping trajectories and net benefits differ for the individual agents that operate the pumping wells identified in Figure 1. The differences are summarized in Figure 5, which shows, for the cooperative and uncooperative strategies, all of the individual agent high-demand shallow aquifer pumping histories, placed next to each other on a single plot with a repeated time axis. These histories show that pumping rates vary slightly across agents for a  given strategy, reflecting differences in the location of their wells. The cooperative pumping rates are lower than the uncooperative rate for all agents at all times.
The western aquifer outflow in Case1 is not constrained and approaches zero near the end of the high-demand 50year simulations, when nearly all the aquifer inflow is captured. If economic conditions and the western boundary head remain constant over time, the western boundary flux will eventually change sign, flowing into (rather than out of) the aquifer. This could have adverse environmental impacts on the downstream environment and could also result in water quality problems for the aquifer if the downstream waters are saline. Such adverse impacts can be mitigated with outflow constraints, such as those discussed in Case 3 below.
The effect of cooperation in Case 1 is to increase net discounted revenue while decreasing pumping, capture, and storage depletion. These simultaneous improvements in both economic and environmental performance are achieved by optimizing pumping over time and across wells to achieve a better balance between demand for water and pumping cost. However, it is important to note that performance metrics such as net benefit, storage depletion, and aquifer outflow depend more on the demand for water and regulatory limits than on the degree of agent cooperation or foresight.
Many analyses of sustainable groundwater management have suggested that, at some point, the cost of pumping groundwater from greater depths will ultimately restrain pumping enough for conditions to reach a new dynamic equilibrium, with stable storage and pumping equal to capture (Negri, 1989). We see a gradual decrease in pumping toward an eventual equilibrium in some of the alternatives plotted in Figure 3. However, when pumping costs are subsidized and/or irrigation water is used to grow high-value crops, it may take a long time to reach an economically constrained dynamic equilibrium. Aquifer outflows and storage can drop substantially in the meantime. This is the case in some of the highly stressed aquifers cited in groundwater depletion assessments (Aeschbach-Hertig & Gleeson, 2012; Bierkens & Wada, 2019).

Case 2. Pumping constrained by well yield limitations
In some situations, well-pumping rates are limited more by well yield than by the cost of pumping (Foster et al., 2015). This effect is addressed in our second case study by considering conditions that activate the well yield constraints described above and in Text S4 in Supporting Information S1. Figure 6 shows how well yield limitations influence the aquifer mass balance for cooperative and uncooperative pumping strategies with high demand in a shallow aquifer. Note that the hydraulic conductivity and specific yield in Case 2 are lower than in Case 1, reflecting a finer aquifer material that gives lower well yields.  In Case 2, pumping activates the well yield constraints at the individual wells after several years into the operating period, resulting in a more rapid decrease in the well pumping rates. This is more apparent in the uncooperative case, where current pumping decisions do not consider future well yield limitations. The effect on saturated depth at the borehole radius at two wells is shown in Figure 7 time series. As the water level drops, the transmissivity decreases substantially, lifts and pumping costs increase, and the unconfined flow nonlinearity becomes increasingly important. Figure 8 shows the individual agent cooperative and uncooperative pumping histories (blue), placed next to one another on a single plot as in Figure 5a. The well yields (brown dashed curves) at some wells are initially above the economically optimal pumping rates but eventually drop in response to decreasing regional heads and near-well transmissivities. When the yield at a given well intersects the pumping rate trajectory, it constrains the pumping rate, and the two rates continue to decrease together. Note that the uncooperative strategy hits the well yield constraint sooner and at more wells. This is the cause of the sharper decrease in total pumping rate shown in the uncooperative aquifer mass balance ( Figure 6: right graph).

Case 3. Pumping constrained by a lower bound on aquifer outflow
In the previous two case studies, agents pumping the example aquifer benefit from the capture of recharge that would otherwise flow out the western boundary to downstream users. Here, we consider how limits on the capture of aquifer outflow might be implemented, and how they might affect net revenue from aquifer pumping. Two major complications arise when considering either voluntary or involuntary constraints to limit capture or, equivalently, to maintain outflow above a specified threshold. These are: (a) the delay between a change in pumping at an agent's well and the corresponding change in boundary outflow (a few years in our example) and (b) the difficulty in assigning responsibility for capture to individual agents.
Cooperative agents make decisions based on complete information about the groundwater system, which includes knowledge of lags in the effects of pumping rate changes and knowledge of the sensitivity of aquifer outflow to pumping by individual agents. This enables them to adjust their pumping rates to ensure that outflows always remain above the desired level. In our analysis, the knowledge available to cooperative agents is contained in the response matrix equations derived from the multi-scale flow model and applicable well yield constraints. They can also include specified lower bounds on aquifer outflow. In effect, the optimization procedure can be used to guide regulators or other decision-makers who might require or recommend individual agent pumping rate limits.  The situation is much different for uncooperative agents, who do not account for either the future impacts of current pumping or the actions of other agents when making pumping decisions. Although the pumping limits obtained from a model-based optimization analysis could conceivably be adopted by uncooperative agents, these agents would then effectively become cooperative by choosing to follow a strategy similar to our cooperative case. A discussion of institutional mechanisms for regulating pumping in uncooperative situations falls outside the scope of this study. We focus here on the impacts of outflow constraints for a cooperative-with-foresight management strategy. These impacts define the maximum net benefit that could be obtained by modifying the pumping decisions of either cooperative or uncooperative agents to ensure that aquifer outflows never drop below a specified target. Figure 9 shows the impact of an outflow constraint for two discount rates over a 25-year operating period. The outflow lower bound (q outmin = 15,000 m 3 day −1 for all years) is 75% of the unpumped outflow value. A hydraulic conductivity of 90 m day −1 , which is somewhat higher than in Cases 1 and 2, facilitates capture of recharge water from the outflow and ensures that the outflow constraint will be activated relatively early in the operating period. The mass balance time series plots show that the cooperative strategy deals with the outflow constraint by pumping aggressively in the early years, before the effects are felt at the western boundary and when the present Figure 8. Yield-constrained annual pumping rates for all nine agents, for cooperative (left) and uncooperative (right) strategies, for silt/clay, high demand, and a shallow aquifer. Figure 9. Outflow-constrained aquifer mass balance flux time series for 3% (left) and 6% (right) discount rates for a cooperative pumping strategy, silt/sand, high demand, and a shallow aquifer. Outflow hits the specified constraint at 25 years.
value of the annual net benefit is high. The total pumping rates fall rapidly, faster for the higher discount rate, until they eventually approach a periodic steady-state with an annual pumping volume equal to the sum of annual aquifer recharge and storage loss (water sources) minus annual outflow (water loss). As the minimum outflow approaches the northern boundary recharge value (no capture) the optimal total pumping rate approaches zero for all years (not shown here), reflecting the fact that pumping can only occur if some of the recharge is captured. A strict outflow constraint greatly limits aquifer pumping and net revenue, compared to the unconstrained case. Figure 10a shows how the cumulative fluxes change when the outflow boundary constraint is activated, around 10-15 years into the operating period. From this time onward, the pumped water for the 6% discount case comes primarily from capture and the cumulative storage depletion levels off at a periodic steady-state value. Figure 10b shows the impact on the pumping rates at individual wells. At the 6% discount rate, the three wells closest to the western outflow boundary stop pumping altogether. The outflow constraint reduces both pumping and benefit substantially, effectively transferring benefit to downstream uses.

Discussion and Conclusions
The three case studies presented above provide quantitative insight that can be used to address our original research questions. We start by considering the impacts of economic and physical factors, scale, and cooperation on pumping decisions. The classical economic controls considered in groundwater depletion studies are the demand for irrigation water, which determines gross revenue, and the cost of pumping and delivering that water. If the marginal demand (the slope of the demand portion of the optimization objective) decreases with the volume pumped, the pumping rate that maximizes net revenue decreases gradually over time, as water levels fall and pumping cost increases. In Case 1, which has no yield or outflow constraints, a higher linear demand term (b 1 ) gives higher pumping and net revenue and is the single most important influence on aquifer depletion.
The physical controls on pumping decisions include aquifer geometry and geology, hydraulic properties, and well spacing. The effects of these controls are conveyed in our analysis through the response matrix elements, which relate annual pumping rates and lifts at individual wells. Unconfined aquifer transmissivity can have a particularly important impact on pumping decisions in areas where the saturated depth is small. This is apparent near our western outflow boundary and also near pumping wells. In Cases 2 and 3, relatively small differences in conductivity can activate well yield or outflow constraints that have significant effects on pumping rates and revenue.
The pumping lift at a given well is defined at a very small spatial scale (the borehole radius) and depends most strongly, at least in the short term, on the local pumping rate. However, over the long term, lift is also affected by the decline in regional water levels that accompanies groundwater depletion. This decline is driven by pumping at other wells and by changes in boundary fluxes. The effect is explicitly represented in our multi-scale groundwater Figure 10. (a) Outflow-constrained cumulative fluxes for 3% and 6% discount rates, cooperative pumping strategy, silt/sand, high demand, and a shallow aquifer. (b) Outflow-constrained annual pumping rates for all nine agents for 3% and 6% discount rates, cooperative pumping strategy, silt/sand, high demand, and a shallow aquifer.
model, which constrains the local head solution to match the regional head on the well circle. As the regional head falls, in response to distant pumping and flux changes, so does the local head. The multi-scale nature of the depletion problem is revealed by time series plots such as Figure 7, which show that local drawdowns are most important over the 6-month pumping period while drops in regional head are important over longer periods.
It is worth noting that some previous studies of the effects of agent decisions on groundwater depletion are based on simplified models that neglect either local drawdowns or the effects of distant boundaries. Examples are the classical bucket model (e.g., Gisser & Sanchez, 1980;Negri, 1989;Provencher & Burt, 1993), which ignores the spatially distributed nature of the aquifer, and models based on superimposed infinite-domain Theis well solutions (e.g., Brozović et al., 2010;Madani & Dinar, 2012a, 2012b. Since these options do not account for inflows or outflows across lateral aquifer boundaries, they cannot consider some of the important issues addressed in our research questions, including the impacts of pumping on downstream water users. These deficiencies have been recognized in more recent groundwater management studies (e.g., Foster et al., 2017;Müller et al., 2017), which are based on finite-domain spatially distributed models that are better able to account for multi-scale effects.
Our study includes another influence on pumping decisions that goes beyond classical demand and cost considerations. This is the degree of cooperation reflected in agent objectives and pumping decisions. Cooperative agents with foresight tend to coordinate their pumping strategies, maximize aggregate benefit, and limit current pumping to reduce future cost. Uncooperative agents without foresight act independently and tend to pump more. These differences are described quantitatively by the different formulations we use for our cooperative and uncooperative optimization problems.
Case 1 focuses particularly on connections between cooperation and performance measures such as revenue and depletion, for different demands and fixed aquifer and cost inputs. Cooperative agents receive higher benefit (as measured by the discounted net present value revenue summed over all agents) and deplete less than uncooperative agents for every option considered. The cooperative advantage reflects benefits gained from coordination of well operations and from knowledge of the future impacts of current pumping. However, our Case 1 results show that differences between cooperative and uncooperative pumping rates, revenue, and depletion are relatively small compared to differences caused by changes in demand (see Table 1). This possibly surprising result deserves further discussion.
An uncooperative agent who ignores long-term regional effects will generally pump more than a cooperative agent who is aware of them. This raises costs as well as benefits, relative to cooperative agent. But both agents will pump at high rates if this is economically advantageous, even if they know that storage will decline and costs will increase in the long term. Cooperation, as we define it here, does not necessarily imply sustainable management. In fact, its economic and environmental advantages tend to diminish as the demand for water increases and there is more incentive for a cooperative group of agents to pump at high rates that yield high short-term revenues that benefit everyone in the group.
Pumping lift and cost may not always be the primary factors that limit pumping. In fact, in situations where energy costs are highly subsidized, pumping cost has negligible impact. This brings us to our next research question, which deals with well yield limitations. In an unconfined aquifer the pumping rate that can be sustained without dewatering the well (i.e., without reducing the saturated depth to zero) gradually falls as storage is depleted and regional water levels decline. This imposes an upper limit on the pumping rate that applies even when pumping costs are small. Figures 6 and 8 show how the Case 2 total aquifer pumping rate falls as yield constraints gradually activate at the nine wells. The yield constraint reduces pumping rates (and net revenue) below the economically optimal values obtained from a classical demand versus cost comparison. Although the pumping rate is still falling at the end of the operating period, the aquifer outflow has dropped to essentially zero, indicating that depletion has reduced benefits for both the pumping agents and for downstream users.
The somewhat lower conductivity and specific yield that apply in Case 2 versus Case 1 cause well yield constraints to limit the Case 2 pumping relative to Case 1, where there are no such limits. However, it is likely that Case 1 yields would also become limiting at some point if operations extended sufficiently beyond 50 years. The effect of yield limitation is especially apparent in Figure 6 for the Case 2 uncooperative strategy, where pumping starts out quite high but must drop more quickly than in the cooperative strategy because well yield constraints are activated sooner. In this particular case, the long-term benefits of cooperation are clear.
Our last two research questions address the effects and regulation of pumping on downstream users, as considered in Case 3. A consistent result in all our cases is that pumping is always satisfied by a combination of diverted recharge and storage depletion (see the case study cumulative volume plots, such as Figure 4). This is a tangible manifestation of the well-known argument against the "water budget myth" (Bredehoeft, 2002;Theis, 1940). The myth suggests that, if pumping demand is less than recharge, it is possible to satisfy demand by extracting water only from recharge (capture) and not depleting storage. Although this may be feasible in an aquifer that functions like a simple bucket, it is not possible to eliminate either capture or storage depletion in our spatially distributed aquifer configuration. This has important implications since capture reduces the water available to downstream users (either groundwater or surface water) and storage reduces the water available to future users of the aquifer. The division between these two sources of pumped water is not easy to control.
To consider this topic further, suppose that we wish to reduce or even eliminate the reduction of aquifer outflow going to downstream users without reducing total pumping from the aquifer, in a situation where pumping is less than natural recharge. Water balance concepts suggest that we may be able to do this by taking more water from storage and diverting less from outflow. In Case 3, we investigate this possibility by imposing a lower bound on aquifer outflow and then deriving the economically optimum pumping rate subject to this new constraint. Figure 9 and Table 1 show that the constrained cooperative optimization problem can satisfy a lower outflow bound set at 75% of the unpumped outflow, but it does this by significantly reducing the pumping rate, which decreases the net revenue to about half of its unconstrained value. Also, Figure 10 shows that storage is still depleted. The results indicate there is more pumping in the first few years with higher discount rates, which put more weight on early revenues, but that the pumping rate is lower in the long term.
Our Case 3 results show that the relationship between pumping, depletion, and outflow is determined by economic factors, the flow system configuration, and the well layout. If we tighten the lower bound so that there can be no diversion of aquifer outflows, the aquifer pumping rate and rate of storage depletion eventually fall to zero. There is no way to pump groundwater over the long term in our configuration without diverting outflow and also reducing storage.
Another interesting aspect of this problem relates to the actual implementation of outflow constraints. As mentioned earlier, cooperative agents can foresee the effects of their pumping decisions and adjust their pumping rates to ensure that outflows always remain above a specified lower bound. By contrast, uncooperative agents without foresight have no way to anticipate the impacts of their pumping decisions on future outflows, even if they wanted to observe outflow constraints. This suggests that foresight, on either the part of agents or regulators, is a prerequisite for effective control of aquifer outflows or, for that matter, other sustainability requirements. In addition, we can speculate that some form of cooperation is also needed to provide the motivation and trust needed to implement outflow or other environmental regulations in practice (Ostrom, 1990).
The integrated modeling-optimization approach described here provides a number of valuable capabilities, but it also has limitations that should be mentioned. One of these is a deterministic formulation, which assumes that recharge and various economic inputs, such as water demand and energy cost, are fixed and perfectly known. A stochastic analysis that accounts for uncertainty and temporal variability, including climate change, would generate results that depend on the methods used to describe the sources of uncertainty. One of the possible advantages of increasing complexity by adopting a stochastic analysis is the ability to consider the effects of extremes (especially droughts) on our conclusions. Other limitations of our analysis include reliance on simplified aquifer geometry and assumptions that recharge to groundwater from direct aquifer precipitation and irrigation return flows are negligible. These hydrologic simplifications are introduced primarily to focus attention on generic issues, such as recharge capture, storage depletion, and well yield limitations. They should be revised as needed for more realistic site-specific analyses.
Our cooperative management strategy assumes that cooperation incurs no transaction costs and has access to perfect knowledge of future aquifer conditions. It would be more realistic to assume that the cooperative strategy operates in real time, optimizing pumping year-by-year, as we have assumed in the uncooperative case. Imperfect knowledge about the future could be accommodated by accounting for prediction uncertainty in the cooperative optimization procedure (Sahu & McLaughlin, 2018b). This real-time version of a cooperative strategy could provide a better assessment of the practical advantages of regulation, as well as cooperation, when environmental and economic conditions are uncertain.
The analysis presented here raises the question of how we should trade off benefits and costs to current aquifer users, future aquifer users, and downstream water users. Storage depletion and outflow reduction are both necessary consequences of aquifer development in our example, even with cooperative management. Although regulatory limits may reduce these adverse impacts, they can also reduce economic benefits and crop production generated from groundwater pumping. It seems reasonable to address this issue by viewing aquifer management from a larger perspective that considers downstream economic and environmental demands as well as current and future uses of the groundwater stock (Gleeson & Richter, 2018). A broader assessment of the depletion problem may reveal opportunities for improved resource use, such as water trading mechanisms, that could redistribute water and associated benefits to everyone's advantage (Wheeler & Garrick, 2020). Analytical methods that combine multi-scale modeling and optimization can play an important role in such assessments.

Data Availability Statement
The inputs, model, and optimization procedure used to generate the results presented in this study are described in detail in the Supporting Information. The software used to generate the results is available from the corresponding author.