Sustainable Groundwater Management in a Two-Cell Aquifer Model

The design of optimal water policies between farmers, municipalities, and groundwater-dependent ecosystems is analyzed in a hydro-economic model with physical interactions between a confined aquifer and a shallow aquifer having a natural drainage. Based on the Pecos Basin case study, we analyze the optimal trajectories of the water tables and water allocation between users and environment flows for the ecosystems. We also explore the consequences for the ecosystem when the water agency uses a one-cell model thinking that it is a correct approximation of the two-cell model. Our results show the importance to consider hydraulic conductivities for the preservation of groundwater-dependent ecosystems.


Introduction
Aquifers constitute a large stock of freshwater. They are critical for human water supply but also for groundwater drainage to their natural outlets, such as rivers, lakes, wetlands, and dependent ecosystems [1][2][3][4]. While climatic change may threaten aquifer recharge [5,6], increased domestic water use with urbanization, the expansion of irrigated agriculture, or the new industrial demand for energy production are actually leading to higher rate of water pumping [7]. For aquifers in interaction with such surface-water bodies, groundwater extraction may entail a decrease in natural drainage and then threaten environmental flows [8]. Determining the conditions of water resources allocation between economic sectors and environmental flows for the dependent ecosystems is a main challenge to ensure a sustainable management of groundwater [9]. Our definition of groundwater sustainability consists in maintaining the aquifer for an infinite time without causing unacceptable hydrologic, environmental, and socio-economic consequences [10]. In particular, groundwater pumping has to conserve discharge for the environment needs [11]. This paper aims at defining what could be such a sustainable management in a multi-use and multi-cell hydro-economic model.
Groundwater management models have been considered under different kind of externalities [12][13][14][15]. The design of water policies has been analyzed for the agricultural sector with one crop [16][17][18] or two crops [19], for the urban sector [20] or for both irrigated agriculture and urban sector [21][22][23]. Part of this economic literature has been analyzed using the "bathtub" or "milk carton" model consisting in a single-cell aquifer with natural recharge as input and natural drainage and pumping as outputs [18]. To deal with conflicting water extraction between two types of users and the preservation of the groundwaterdependent ecosystems, a more complex hydro-economic model has to be considered. For that purpose, we develop a two-cell aquifer model with a confined aquifer and a shallow aquifer which has a natural drainage to a dependent ecosystem. Such a configuration refers to the initial modelling of the Roswell groundwater basin and the Pecos River by [24] for which data are available. It is worth to mention that most of the economic literature after [18] used the one-cell aquifer version, often without natural drainage, under the label of Pecos Basin and not the original model. This is not without raising many questions concerning the hydraulic coherence of the system and the consequences on the design of water policies. We thus develop a more general than [24] and consider a stylized multi-use and multi-cell hydro-economic model with two groups of users, municipalities and farmers, having their own demand function. Urban water use is extracted from the confined aquifer while water dedicated to agriculture is pumped from the shallow aquifer. Modelling such interactions is important for the design of water pumping policies which consider that groundwater and surface water are not separate resources and are functionally interdependent by their hydraulic connection.
The hydrological literature on the conjunctive use of surface and groundwater at a basin scale exists but is scarce. Some economic models have considered several connected bathtub aquifers, several sources of water between surface water, rainwater, and aquifer water [25][26][27][28][29]. These models often focus on irrigated agriculture with farmers maximizing their payoffs. In terms of hydro-economic modelling [24], consider a confined aquifer which is located below an unconfined aquifer with vertical flows governed by Darcy's law implying that water flows upward due to the water pressure. They also assume that the confined aquifer has a unconfined part with a natural recharge distinct from the shallow or unconfined aquifer. It differs from the literature on multi-cell models assuming water flow exchanges between lateral and not vertical aquifers as in [30][31][32]. [30] consider a model of several interconnected cells but water flows is limited to between adjacent cells with one user on each cell while [31,32] consider two adjacent cells. Even with lateral connections, water flows between the two cells are proportional to the difference in stock levels and thus are governed the Darcy's law. By contrast, [33] also consider a system of confined and unconfined aquifer but assume perfect transmissivity between the two aquifers. [34] isolate the water pressure externality in the particular case of an artesian aquifer which is a confined aquifer that recharges from a more elevated unconfined aquifer without any pumping cost externalities. They also assume perfect transmissivity between the recharge area and the confined aquifer. However, the modelling approach of [24] do not take into account spatial externalities as in [35] which consider that pumping can induce cones of depression around individual wells based on Theis's argument [36]. It implies that when users or wells are closed together, the overlap of depression cones may reinforce the pumping externality. [37] considers the issue of saline intrusion in a hydrological system with spatial externalities in which the confined aquifer recharges water from the unconfined aquifer and discharges into the sea. [37] shows that the seawater externality may exceed the pumping cost externality near the coast. [38,39] investigated more complex models with real-world cases but they have not been extended in the economic literature.
This article aims at developing an intermediate stylized model of a synthetic basin dealing with economic, hydraulic, and ecological interactions between users, between aquifers, and between aquifer and a river which determine the amount of environmental flows for the ecosystem. When natural drainage sustains environmental flows for a river, the decrease in the water table reduces the natural drainage and creates damages to the ecosystems [14,40]. It becomes critical for ecological safety that the water table of the shallow aquifer remains above a critical level to ensure that water flows from aquifers to the river [41][42][43]. When the flow is in the opposite direction, the aquifer is fed by the river and it may have disastrous environmental consequences in terms of biodiversity losses for instance. The existence of these environmental flows is of importance for a water agency aiming at maximizing a welfare function, defined as the sum of the net benefits of all the agents plus the ecosystem damages, and then internalizing both the pumping cost and the environmental externalities.
Our contribution to the literature in terms of groundwater management concerns the choice for a water agency of using the one-cell instead of the two-cells model to determine his optimal policies. If conditions ensuring that the two models are not equivalent from a hydrological point of view, then it is important to analyze the consequences for the water agency to use, as an approximation of the initial two-cell model, a one-cell model as if it was the correct model. In this case, it is important to see whether the water agency can maintain positive environmental flows and preserve the ecosystem. We know that such an objective requires to determine how much the society has to value in monetary terms the ecosystem services, sustained by groundwater flows, to preserve them [44][45][46]. Hence, the value of the marginal ecosystem damage which determines the weight of the environmental externalities may differ according the hydraulic models, meaning that the water agency has to allocate differently water pumping between farmers, municipalities, and the ecosystem. Section 2 presents the aquifer model by distinguishing the one-cell model with the two-cell model and the equivalence conditions between them. Section 3 analyzes in a dynamic game how a water agency maximizes a welfare function including ecosystem damages due to aquifer depletion. Section 4 illustrates our results based on the hydraulic system of the Pecos Basin. Section 5 concludes. All the proofs are given in the Appendix.

Groundwater Model
This section presents the hydraulic components of the model and in particular the equivalence conditions between the one-cell and two-cell models are analyzed.

From One-Cell To Two-Cell Aquifers
Time is indexed by t ∈ ℝ + . The continuous time hydrological dynamics is presented for both one-cell and two-cell models. When a two-cell hydrological model is considered, the dynamics of an aquifer k ( k = 1, 2) describes changes in the water table, measured by h k ∈ [0, h max k ] where h max k stands for the maximum level of the water table. We also introduced a minimum value for the water table h min k due to the presence of natural drainage W f consisting in environmental water flows for river, lakes, or groundwater-dependent ecosystems. When only one aquifer is considered (Fig. 1), subscript k is omitted. When two aquifers are considered (Fig. 2), h 1 is the water table of the confined aquifer and h 2 the water table of the shallow aquifer. Water tables increase with the constant natural recharge R k > 0 and decrease with water pumping and the natural drainage. The hydraulic connection between the two aquifers is given by the flow D 1→2 . As explained in the introduction, these one-cell and two-cell models are based on the hydraulic features of the Pecos Basin by [18,24]. We only depart from these models by considering water extraction W a and V m dedicated to irrigated agriculture for n a identical farmers and to water consumption for n m identical municipalities. Then, W a = ∑ n a i=1 w i where w i stands for the individual water extraction of farmer i ∈ [1, .., n a ] and V m = ∑ n m j=1 v j where v j stands for the individual water extraction of municipalities j ∈ [1, .., n m ] . A proportion i with i ∈ {a, m} of the water is assumed to return to the aquifer or in each aquifer where 0 < i < 1 stands for the return flow coefficient. A k is the area of aquifer k and S k the storativity coefficient.
Based on Fig. 1, the dynamics of the one-cell aquifer is where ̇h stands for d dt (h(t)) . According to [47], the natural drainage describing the stream-aquifer flow is based on a linear conductance model. This assumption is mainly accepted in the literature [48,49]. We have We assume that in the pristine initial conditions before pumping starts, natural discharge is equal to the natural recharge R = W f . Under this condition, the conductance and constant coefficient k f can be rewritten as follows (2) shows that for h ≥ h min water flows from the aquifer to the river but for h < h min the aquifer is fed by the river. Such a case may have disastrous environmental consequences in terms of biodiversity losses and is in contradiction with our definition of sustainability. This critical water table h min can be interpreted as a tipping point in this hydrological system. We assume that in this configuration the river does not become completely dry, implying no regime shift (see [50,51] for examples of regime shifts).
In the one-cell model, water extractions for farmers and municipalities are perfectly substitutes and impact symmetrically the natural drainage.
From Eqs. (1) and (2), the dynamics can be rewritten as with the following notations (1-µ m )·V m Fig. 1 One-cell aquifer model with two uses, one recharge and a natural drainage In the two-cell model, we assume that the municipalities extract water from the confined aquifer h 1 while extraction dedicated to irrigated agriculture is from the shallow aquifer h 2 . We also consider that the natural drainage, now denoted by W n , is from the shallow aquifer. It implies that water extraction by farmers in the shallow aquifer impacts directly the natural drainage while the impact of water extraction in the confined aquifer is indirect through the hydraulic connection between the two cells. E measures the evaporation by phreatophytes in the dynamics of the shallow aquifer 1 . The dynamics is given by The flow D 1→2 between the two aquifers obeys Darcy's Law where k 12 > 0 is the leakage proportionality coefficient The flow W n from the shallow aquifer to the river is determined by As in the one-cell model, we assume that in the pristine initial conditions before pumping starts, natural discharge is equal to the natural recharge

Equivalence Between the One-Cell and Two-Cell Aquifers
Due to its larger size and the interactions involved, the twocell model is more cumbersome to analyze than the one-cell model. Taking into account the complexity of the two-cell model will therefore only be of interest if the trajectories of the two models are significantly different. For that purpose, we aim at looking to which extend the one-cell model can be a good approximation of the two-cell one, meaning that, for given water extractions W a and V m , trajectories of the two systems are very close. Three distinct situations in the two-cell model have been identified depending on the values of specific coefficients which ensure that the two-cell model behave as a one-cell one. Results are presented in Propositions 1-3. Proofs given in the Appendix are based on properties of perturbation theory [52,53]. The three identified cases are (i) the natural drainage can be considered as negligible, (ii) the storage capacity of one aquifer (i.e., the product of the area of the aquifer and the storativity coefficient, AS) is small with respect to the other aquifer, and (iii) leakage rate k 12 is large.
In Proposition 1, we consider natural drainage parameter is negligible. To do so, we introduce a parameter on the natural drainage parameter. We then analyze what happens when coefficient tends towards zero. The model appears to be a regular perturbation of a simpler dynamical system. When the natural drainage is small, it means that whatever the level of groundwater table, there are few impact on the ecosystem in both the one-cell and two-cell models. Proposition 1 then proves that in that case, trajectories of the two-cell model are indistinguishable from the ones obtained in a one-cell model, in which the water table height is the weighted mean of the water table in each cell and where the recharge is the total recharge that impacts the whole groundwater. Formally, Proposition 1 is given as follows.
Proposition 1 Assume k n = n with 0 < ≪ 1 , which means that the natural drainage is negligible. Let h be defined as h = Then, the dynamics of the mean water table h obtained from the two-cell dynamics given by system (5) with coefficients y a and y m defined in Eq. (4) is It is a regular perturbation, for which solutions are close to the reduced dynamical system which corresponds to the dynamics of a one-cell aquifer.
Equation (11) is called the dynamics of the reduced model obtained when tends toward zero. Figure 3 illustrates the discrepancy between the two dynamics obtained for the reduced model and for the two-cell model in terms of water table and in volume, when the difference of the

3
water table is multiplied by the water storage capacity of the aquifer. We consider two values for parameter k n with k n = 0.1 and k n = 2.8783 which corresponds to the true value for the Pecos Basin (see Tables 1 and 2 in Subsection 4.1). We also consider constant extraction water with V m = 268, 64 Mm 3 and W a = 262, 89 Mm 3 . Figure 3a shows that the smaller k n , the closest the trajectories of the onecell model to those of the two-cell model. Moreover, the discrepancy increases with k n . For k n = 2.8783 , the difference is about 0.5 m in terms of water table but at the scale of the aquifer in volume (multiplied by the term AS), the discrepancy becomes important as shown in Fig. 3b. Now singular perturbation tools are considered by introducing a coefficient on the water table of the confined aquifer which is a state variable of the system. We assume the storage capacity of the confined aquifer, A 1 S 1 , is very small compared to the shallow one, A 2 S 2 . In that case, the dynamics of the two-cell model can be rewritten as a twotime scale dynamics. Using Tikhonov's theorem [52], it can then be shown that the trajectory quickly converges to a manifold (called the slow manifold and denoted L), and then the dynamics on the slow manifold is the one of a one-cell model. It gives Proposition 2: Then, the dynamics of the two-cell aquifer is given as follows.

For t large enough, the water table of the confined aquifer satisfies
where the dynamics of the non-confined aquifer can be approximated by the reduced model which is a one-cell aquifer dynamics, namely: Corollary 1 gives the formal formulation of Proposition 2.

Corollary 1 Let us assume that
Then, the dynamics of the two-cell aquifer is a two-time scale dynamics where t is the slow-time. Slow time t dynamics is given by system (12).
At fast time defined as t = , the dynamics quickly converges in the neighborhood of the attractive slow manifold, defined as  where the slow motion takes place satisfying Eq. (14). Figure 4a illustrates the phase diagram of the singular dynamics of the two equations of system (12). The phase diagram shows that the faster the convergence to the slow manifold, given by the solid line (in red) and by Eq. (13), the lower A 1 S 1 is or equivalently . It also can be shown that the convergence to the slow manifold is very quick. Figure 4b shows the trajectory of the water table of the shallow aquifer h 2 for different values of A 1 S 1 of the confined aquifer. As proved in Proposition 2, the smaller A 1 S 1 , the closer the trajectories to the one obtained in the reduced model given by Eq. (14). We consider two different values of A 1 S 1 with A 1 S 1 = 25 and A 1 S 1 = 101 which is the true value of our case study the Pecos Basin. While A 1 S 1 = 101 corresponds to 22.6% of the value A 2 S 2 for the non-confined aquifer, the lower value A 1 S 1 = 25 is equal to 5.6% of the value A 2 S 2 for the non-confined aquifer. Figure 4b shows that for A 1 S 1 = 101 the one-cell dynamics cannot be used as an approximation of the two-cell model given by Eq. (12). However for a lower value of A 1 S 1 in terms of A 2 S 2 , the paths of the water tables of the one-cell aquifer and the two-cell aquifer are very close over the simulation period.
Here again, a singular perturbation theory is used to study the dynamics of the two-cell aquifer when we consider the leakage coefficient k 12 between the two aquifers.
Proposition 3 shows how the trajectories of the two-cell model can be approximated by solutions of a one-cell model for large leakage coefficients.
. Then for t large enough, the water tables of the confined and non-confined aquifer are close, namely: and the dynamics of the two-cell aquifer is close to the dynamics of the following one-cell aquifer The formal expression of the dynamics using Tikhonov's theorem is expressed in Corollary 2.
. Then, System (5) can be rewritten as the following slow-fast dynamical system.
Let us consider initial conditions h 0 , h 2 (0) . At fast time defined as t = , the dynamics quickly converges in the neighborhood of the attractive slow manifold, defined as for the reduced model (Eq. 14) in solid (red) and for system (12) for

3
where the slow motion takes place satisfies Equation 17 is called the equation of the reduced model. Fig. 5a illustrates this singular dynamics for two different values of the leakage coefficient, k 12 = 33.395 which is the true value of the Pecos Basin and a new larger value k 12 = 100 . The phase diagram shows that the larger k 12 , the faster the convergence to the slow manifold (solid line (red) of Fig. 5a). Fig. 5b shows that the larger k 12 , the closer the trajectories are to the ones of the one-cell model.

Groundwater Management
We assume that there exists a water agency aiming at maximizing a welfare function defined as the sum of the benefits of the farmers and municipalities and an ecosystem damage due to aquifer depletion under the constraint given by the aquifer dynamics. We first describe the behavior of farmers and municipalities and the environmental externalities and derive the social optimum achieved by the water agency.

Farmers and Municipalities
As in [18], water demand functions are a linear function of the price p w and differ between the farmers and municipalities. W a being the total demand for farmers and V m being the water demand for municipalities, we have The benefit of a farmer i becomes In a similar way, the benefit of a municipality j is The marginal cost of pumping depends on the level of the water table such that for i ∈ {a, m} It shows that the marginal cost of pumping increases with the intercept of the pumping cost function c i 0 which gives the maximum average cost of extraction and decreases with the slope of the pumping cost function c 1 which depends on the water table. We can define h max k as the maximum water table elevation associated with the natural hydrologic equilibrium of the aquifer that occurs when groundwater reserves reach the storage capacity of the aquifer. Hence for that value of the water table h k = h max k , the marginal cost of pumping is zero. Then from Eq. for the reduced model (Eq. 15) in solid (red) and for system (16) for k 12 = 100 in dot (black) and k 12 = 33.395 in dash (blue) 1 3 [54]. These coefficients can differ between the municipalities and the farmers in the two-cell model since water is extracted from different aquifers. However in the one-cell model, we assume that they are equal c i 1 = c 1 since we assume that both users face the same hydrogeological conditions. In that case, cost heterogeneity can be introduced in the fixed cost due to different delivery charges [20].
In the one-cell model, the benefits derived for a farmer i from pumping a water quantity w i and a municipality j from pumping an amount v j are given by with We obtain similar expressions in the two-cell model. The payoff functions for a farmer i pumping an amount w i in the shallow aquifer and for a municipality j pumping an amount v j in the confined aquifer are with the new coefficients

The Ecosystem Damages
According to [40], the ecosystem damages function due to the depletion of the aquifer is assumed to take the following form The first term of Eq. (25) consists in the cost of the aquifer depletion measured by the capture variable defined by [55] as the decrease in drainage plus the increase in recharge (not considered here since R is constant). This is of interest for damage to ecosystems associated with consumptive uses (river base flow, transpiration by phreatophytes). Coefficient > 0 measures the cost of damages to the ecosystem for each cubic meter of depletion. The second term of Eq. 25 refers to another kind of damage cost related to the difference between the maximum (initial) and the current water levels as in [56]. It is mostly relevant for nonconsumptive uses of groundwater (for example, to avoid subsidence). Coefficient > 0 is also a measure of the cost of damages to the ecosystem for each meter of depletion (or height of water level drop). Using Eqs. (2) and (25) can be rewritten as follows: In the two-cell aquifer, the ecosystem damage is given by which can be rewritten as

The Social Optimum
In the configuration of one aquifer and two uses, the water agency aims at maximizing the sum of the benefits of the n a identical farmers and the n m identical municipalities and also taking into account the environmental damage given by Eq. (26) 2 . The program of the water agency can be written as To ensure an interior solution, we introduce the following assumptions Assumption A 1a ∶ r − (blaya+ l m y m )( +e) Assumption A 1b ∶ e + c 1( l a y a +l m y m )( +e) +e+c 1( l a y a +l m y m) − c 2

Proposition 4 Under Assumptions A 1a and A 1b , the water table dynamics for the social optimum is given by
where h SO is given by and the eigenvalue with negative real part given by > 0, 2 The social optimum can be implemented with market-based instruments as water market between municipalities and farmers in two connected aquifers. Such market-based instruments have already been analyzed between two groups of farmers [19] or between one group of heterogenous farmers [57][58][59].

3
It can be shown that the water table at the steady state does not depend on the number of users.
As a particular case, assumed that there exists only one user, no natural drainage (implying e = 0 ) and also no environmental damage d 1 = 0 , h SO can be rewritten as: which is the same expression as the one obtained by [51,54]. An increase in water demands reduces the water table while this is the opposite for higher fixed costs. An increase in the slope of the damage function implies a higher water table. The assumption y a = y m = y implies that the return flow coefficients are the same for the two uses ( a = m = ). It can be shown that under this assumption, the above assumptions A 1a and A 1b are fulfilled.  In the two-cell aquifer configuration, the program of the water agency becomes under the dynamics with the initial conditions h The resolution of the water agency program gives the following proposition: The next section illustrates the dynamics obtained in the social optimum for the one-cell and the two-cell aquifer models.

Numerical Illustrations and Management Policies
This section presents the calibration of the one-cell and twocell models, then the impacts of the different coefficients on the key variables of the model at the steady state and the management policies to be implemented to preserve the ecosystems. The management scenarios are numerically tested using data on the Pecos Basin. Data from [24] for the two-cell model and [18,47] for the one-cell model are rather old but they are the only available ones. Concerning the economic calibration of the parameters, we follow [29] which have converted the 1970s dollar into the equivalent 2012 value (multiplied by 5.49), and then into euros, with an exchange rate of 0.778€/$.

Parameter Values
Tables 1 and 2 provide the hydrogeological data used for the two-cell model and the one-cell model.
As the values of h max and h min are the same as the confined aquifer, it implies that the conductance coefficients k f and k n are the same. As explained in Subsection 2.2, conditions ensuring the equivalence between the two models are not satisfied. The natural drainage is not negligible, the storage capacity of the two aquifers are significant, and the leakage coefficient between the two cells is not large. It means that the hydrologic properties of the two models are different. As pointed out in the introduction, this difference is important if we consider that the manager wrongly thinks that the one-cell model he uses is a correct approximation of the two-cell model. It should be noted that when [18] use their representation of the two-cell model by a one-cell model, they also assume that the natural drainage was zero, implying the equivalence between the two models.
The Pecos Basin case study involves only one group of farmers pumping both in the shallow and the confined aquifer at different costs. As we introduce two groups of users, farmers and municipalities, we assume that each group has a specific water demand function and can only pump in one sub-aquifer (see Table 3). However to keep the same water extraction pressure, we assume that the sum of the two water demands is equal to the aggregate demand in [18], i.e., W = g − kp w with W = W a + V m , g = g a + g m , and k = l a + l m . To simplify, we consider that the two water demand functions are the same but costs are different as in [24]. It can be checked that  [60,61]. In a different case study, [62] have considered for the coefficient a value of 0.5€ per m 3 and for a value of 100000€ in [56]. However, such data are not available for this case study and another approach must be followed as shown in the next Subsection 4.3.

Steady-state Comparative Results
As shown in Corollary 3, the water table and the amount of water extraction by the users do not depend on the number of users. It means that when the number of players varies, the global amount of water extraction for each group remains the same but more users inside a group mean less water for everybody. In the one-cell aquifer and at the steady state, Table 4 shows the impacts of a shift on the coefficients on the variable of interest. An increase of the municipal water demand g m shows that more water is allocated to municipalities and less to farmers. However, the net impact on the aggregate extraction is positive and implies a lower water table. A rise in the slope of the water demand function l m or in the fixed cost c m 0 implies an increase in the water table due to a lower water extraction from municipalities. It can be shown that even if water extraction for the farmers increases, the overall impact on the total of extracted water remains negative. As we assume that both groups of users have the same cost coefficient c 1 , the impact of a rise in c 1 decreases the marginal pumping cost implying higher water extraction for both users at the expense of a fall in the water table. Lastly, a rise in d 1 means less water for users and a higher water table.
Results also show that the impact of a rise in the natural recharge R depends on the presence of the environmental externality in the welfare function of the water agency. We refer to notation SO for the social optimum when the pumping cost externality is internalized but not the environmental externality ( d 1 = 0 ) and notation SOE when the environmental externality is also taken into account ( d 1 > 0 ). Table 5 shows that a rise in the natural recharge R always increases the water table in the two configurations but while water pumping rises in the SO case, it decreases in the SOE case. It shows the impact of the environmental externality since a rise of the natural recharge also means a higher natural drainage through the higher conductance coefficient k f and it also increases the weight of the environmental externality given by d 1 in the welfare function of the water agency.
In the two-cell aquifer, Table 6 gives the changes on the main parameters on the variable of interest. A positive shift on the municipalities' water demand ( g m ) implies a transfer of water from farmers to municipalities. As total water extraction increases, the water tables of the two aquifers decrease. An increase in the slope of the municipalities' water demand ( l m ) or a rise in the fixed pumping cost implies a transfer of water from municipalities to farmers. The net impact on water extraction is negative and implies a rise of the water tables. When the slope of the pumping cost function ( c m 1 ) increases, it reduces the marginal cost of pumping and increases water extraction for municipalities. The decrease in farmers extraction is not enough to offset the municipalities' increase, leading to a higher total extraction and a fall of the water tables.
A full analysis of these comparative results after a 10% shift in several parameters is provided in Tables 12 and 13 in Appendix 7.8.

Water Management Policies for the Preservation of the Ecosystem
As conditions ensuring that the two models are not equivalent from a hydrological point of view, we now assume that the water agency uses the one-cell model as if it was the correct model and not the two-cell one to determine his optimal management strategy. In particular, the water agency aims at ensuring positive environmental flows and preserving the ecosystem. We know that such an objective requires to determine how much the society has to value in monetary terms the ecosystem services sustained by the groundwater. In our Table 5 Impact of a rise of the natural recharge R in the onecell model    (26) and (27) and associated coefficients and . Due to the lack of data, we first have determined the values of theses coefficients in such a way that the steadystate water table in the one-cell aquifer remains above the critical threshold h ≥ h min (see Table 7). Table 8 gives the steady-state values in the SO configuration (when d 1 = 0 ) and SOE one (when d 1 > 0 ) where Q  Then, based on these value for coefficients and , our objective is to see if the condition h 2 ≥ h min 2 is satisfied in the two-cell aquifer model. Table 9 shows that the water table h 2 at the steady state remains below the critical value h min 2 whether we assume or not cost heterogeneity between farmers and municipalities. Case (a) refers to identical cost coefficients ( c a 0 = c m 0 and c a 1 = c m 1 ) while case (b) refers to heterogenous cost parameters ( c a 0 ≠ c m 0 and c a 1 ≠ c m 1 ) The trajectories of the water table and the pumping for both uses in the two configurations, SO and SOE, are given in Fig. 6 in the one-cell model and in Fig. 7 in the twocell model. We also assume that the initial conditions for the water tables are their maximum values. In both models, the number of users in each group is identical and given by n a = n m = 100.
Results show the differences in the trajectories of the water table in the one-cell and two-cell models (both the confined and the shallow aquifers). The environmental externality is measured by the difference between SO and SOE. As explained above, the values of the costs of capture and ecosystem damages have been chosen to maintain the steady-state water table in the one-cell model at its critical value h min in the SOE configuration. A higher weight of the environmental externality implies a higher water table and a lower amount of water extraction. Our results show that such values are not enough to ensure this condition in the two-cell model. With the same values of the costs of capture and damages, the water table of the shallow aquifer at the steady state is lower than its critical value h 2 < h min 2 meaning that the river is fed by the shallow aquifer in the SOE case. This result holds when the cost functions of the two groups of users are identical and adding cost heterogeneity increases the difference between h 2 and h min 2 , from 1.3 to 2.2m at the steady state. Such a difference relies on the hydraulic conductivities which impact the design of groundwater management. In the one-cell model, farmers and municipalities' water extractions impact symmetrically the water table of the aquifer while in the two-cell model the impact of water extraction of municipalities on the shallow aquifer is indirect. It results a higher water extraction for municipalities ( 1.523Mm 3 and 1.557Mm 3 instead of 1.458Mm 3 ) and a slightly lower water pumping for irrigation ( 1.445Mm 3 and 1.446Mm 3 instead of 1.458Mm 3 ) in the long term. However during the 300 simulation period, Fig. 8 shows that water   extraction in the SOE case for both uses is higher in the twocell model compared to the one-cell model.
In the two-cell model, it suggests that a higher monetary value is needed to ensure the preservation of ecosystem services sustained by environmental flows. It implies that if the water agency uses the one-cell model, which is a wrong simplification of the two-cell model, the value of environmental externality will be underestimated in the design of optimal policies in the two-cell model, implying too much water pumping for farmers and municipalities. To achieve the objective of positive environmental flows h 2 ≥ h min 2 in the two-cell aquifer, the values of coefficients and have to be respectively higher as shown in the next Table 10. Now we look at the limit case when k 12 becomes large ( k 12 → ∞ ). As shown in Subsection 2.2, such a condition ensures that the one-cell model behaves like the two-cell model. Results in Table 11 show that for the set of environmental damage coefficients defined in Eq. (7) the water table tends to its critical value h min in the SOE case.
Our results show the importance of taking into account the true hydrogeological model, in this case the two-cell model, in the design of environmental policies. When the one-cell model is a wrong approximation of the true model, total water pumping appears too excessive to preserve positive environmental flows.

Conclusion
The paper aims at studying how the introduction of hydraulic conductivities in a stylized two-cell hydro-economic model impacts groundwater sustainability when different groups of users, farmers, municipalities, and ecosystems are in competition for the resource. The social optimum water policies have been analyzed when environmental externalities are internalized by a water agency.
Our result shows that the two-cell model as initially described by [24] could not be approximated by a one-cell model due to the hydrogeological system of the Pecos Basin. This is of particular importance when an environmental externality is added implying that the aquifer is open with natural drainage [40]. It should be noted that when [18] use the onecell model, they assume a closed aquifer with no natural drainage. Results show that if the water agency uses the one-cell model as if it was a correct representation of the initial two-cell model for the design of his water policies, his objectives in terms of preservation of ecosystem will not be achieved. In particular, we show that the value of the weight of the environmental externality which ensures the preservation of positive environmental flows in the one-cell model is too low compared to what was required for the two-cell model. It means that the water agency allocated too much water for farmers and municipalities at the expense of the ecosystem. It suggests that the use of the one-cell bathtub aquifer model even when natural drainage is explicitly taking into account may appear too simple to deal with groundwater management and sustainability. It justifies the need to consider hydraulic conductivities and multi-cell models.
Overall, our methodology can be used to handle more complex systems and different hydraulic configurations. Potential extensions include the analysis of switching regimes in the dynamics of the water table. We have explicitly assumed that when the water table of the shallow aquifer was below a critical level, water flows from the river to the aquifer. However, the case of a dry river was excluded. Relaxing this assumption will be addressed in future research endeavors. Author Contributions Emmanuelle Augeraud-Veron has been involved in all the stages of the paper (mathematical modelling, simulation, redaction). Jean-Christophe Pereau has been involved in all the stages of the paper (mathematical modelling, simulation, redaction).

Funding Information
The paper has the support of CNRS-CSIC International Research Program 'ALLIES' 2020-2023 GREThA-IPP.

Availability of Data and Material Not applicable.
Code Availability Not applicable.

Conflict of Interest
The authors declare no competing interests.

3
with y a , y m defined in Eq. (4). Dynamics is a regular perturbation of Eq. (40). Thus, the trajectories of System (11) and the dynamics of the two-cell model are very close when k n is small.

Proof of Proposition 2
Let us consider the following two-time scale dynamics written in slow time t where notation ̇x stands for d dt (x(t)).
with initial condition (h 0 1 , h 0 2 ) . System (41) has two limit systems. At fast time , defined as t = , system (41) is a regular perturbation of the unperturbed system At fast time , h 2 is constant, equal to h 0 2 , and h 1 varies quickly and is approximated by the solution of the boundary layer The slow manifold is defined as is attractive and thus the theorem of Tikhonov [52] can be applied. Let us denote Letting (h 0 1 , h 0 2 ) be an initial condition, as the slow manifold is attractive, according to Tikhonov's theorem, the solution of system (41) quickly converges to the neighborhood of the slow manifold according to Tikhonov's theorem. And then, on the slow manifold, the slow motion takes place. The slow motion is given by Substituting H 1 (h 2 ) gives Eq. (14). Then, Proposition 2 holds.

Proof of Proposition 3
Let us consider the following slow-fast dynamical system with initial condition h 0 , h 2 (0) .
At fast time defined as t = , the dynamical system is equivalent to the unperturbed dynamics Slow manifold is thus defined as As the slow manifold is attractive, Tikhonov's theorem applies, and thus a trajectory with initial condition h 0 , h 2 (0) quickly converges in the neighborhood of the slow manifold, and the dynamics on the slow manifold, where h 2 = h , the dynamics is then given by Eq. (15) in the main text.

Proof of Proposition 4
The Hamiltonian relative to the manager optimization problem is given by where stands for the co-state variable. The optimality conditions are together with the transversality condition lim t→∞ e − t h(t) (t) = 0. We first assume that the control constraints are not binding. Substitute the optimal expressions of the water extractions in the dynamics of the aquifer yields At the steady state, we obtain The steady-state water table is positive h SO > 0 when the numerator and denominator are positive, under Assumption A 1a ∶ r + e + c 1 l a y a + l m y m − bl a y a + l m y m ( + e) −c 1 l a l m y a − y m by m − y a > 0 and Assumption A 1b ∶ e + c 1 l a y a + l m y m ( + e) + ec 1 l a y a + l m y m − c 2 1 l a l m y a − y m 2 > 0. The steady-state value of the adjoint variable is given by: The dynamics of linear system (43)  the determinant of the Jacobian matrix which is negative according to Assumption A 1b . The characteristic polynom Δ(X) is given by Δ(X) = X 2 − X + . As < 0, the discriminant is positive; thus, there exist two real roots solution of equation Δ(X) = 0. These eigenvalues are explicitly given by = ± √ 2 − 4 but the one which satisfies the transversality condition is

Proof of Proposition 5
The optimum, at the steady state, gives the amount of water extraction for farmers  with which can be rewritten as and for the municipalities

Proof of Corollary 3
From the expression of h SO given in Eq. (30), it is immediate to see that h SO is independent of n m and n a . Using the expressions of in Eq. (35) and w SO in Eq. (33), w SO is also independent of n m . Using the expression of v SO in Eq. (34), a simple computation shows that v SO n m = −1 . To simplify the expressions, we assume y a = y m = y.
The derivatives of h SO with respect to coefficients b and are Based on the expressions of b and , we get We now consider the impact of c 1 . Using the expression for h SO we have Thus the sign of h SO c 1 is the same as the sign of the following term which is negative as b and < 0.

Proof of Proposition 6
The Hamiltonian of the system is defined as The proof is decomposed in several steps: In step 1, we derive optimality conditions. Step 2 is devoted to the computation of the long run optimal solution ( h SO 1 , h SO 2 ). In step 3, we prove the existence of saddle-path configuration where the stable manifold is two dimensional. In step 4, we write the explicit dynamics of the solution, using Riccati decomposition.
Step 1: Optimality conditions are   and   The optimality conditions together with the water table  dynamics give Step 2: Water tables h  = + e 0 − n m y 1 2 − e 1 + n m with Step 3: The Jacobian matrix derived from the optimality conditions dynamics is Characteristic polynom is where M 2 and M 3 are the sum of all diagonal second-and third-order minor of J. Properties of these kind of Jacobian matrix are well known [63], and it can be easily checked that −M 3 + M 2 + 3 = 0.
Thus   All the terms are positive except the last term. Assuming 4f y 1 y 2 e 0 e 2 − f y 1 e 1 + y 2 e 0 2 > 0 e n s u r e s t h a t det (J) > 0 and as K < 0 , there exists two roots with negative real part, denoted by Step 4: We now turn to the explicit solution of the optimal dynamics The optimality conditions together with the water table dynamics give Letting H � = h 1 , h 2 � and Λ � = [ , ] � , the system of the optimal conditions and the water table dynamics can be rewritten as h 2 = − y 2 a n a 2d + e 1 h 1 + fy a n a 2d − e 2 h 2 + r 2 + by a n a 2d .

Steady-state Comparative Results
Based on the steady-state values in the one-cell and twocell model which correspond to our benchmark case, we analyze the impact of a 10% shift in several parameters concerning the number of players, the water demand function, and the pumping cost function. Table 12 gives the impact of a shift of 10% on the different coefficients at the steady-state variables for the one-cell model. The first number is the absolute change and the second number is in % with respect to the initial situation. For instance, an increase of the municipal water demand g m by 10% (from g m = 290 to g m = 319 ) shows that the water consumption for the municipalities increases by 10.75% (or +27.07Mm 3 ) while the farmers' consumption decreases by 0.661% (or -1.92Mm 3 ) in the SO configuration. Table 13 gives the impact of a shift of 10% on the different coefficients at the steady-state variables for the twocell model. The first number is the absolute change and the second number is in % with respect to the initial situation with cost heterogeneity.