Ecosystem Shifts: Implications for Groundwater Management

Freshwater ecosystems provide a large number of benefits to society. However, extensive human activities threat the viability of these ecosystems, their habitats, and their dynamics and interactions. One of the main risks facing these systems is the overexploitation of groundwater resources that hinders the survival of several freshwater habitats. In this paper, we study optimal groundwater paths when considering freshwater ecosystems. We contribute to existing groundwater literature by including the possibility of regime shifts in freshwater ecosystems into a groundwater management problem. The health of the freshwater habitat, which depends on the groundwater level, presents a switch in its status that occurs when a critical water level (‘tipping point’) is reached. Our results highlight important differences in optimal extraction paths and optimal groundwater levels compared with traditional models. The outcomes suggest that optimal groundwater withdrawals are non-linear and depend on the critical threshold and the ecosystem’s health function. Our results show that the inclusion of regime shifts in water management calls for a reformulation of water policies to incorporate the structure of ecosystems and their interactions with the habitat.


Introduction
A recent report from the United Nations (IPBES Global Assessment Report on Biodiversity and Ecosystem Services) has warned of the unprecedented rapid deterioration of the health of ecosystems worldwide (Díaz et al. 2019). This statement provides evidence that environmental regulations for protecting natural resources and the defense of ecosystem's health and functioning are still far from being accomplished. Human activities are the main drivers of ecosystem damage, causing several problems in the functioning and health of ecosystem organisms, loss of ecosystem resilience, loss in biodiversity, and habitat destruction (Lu et al. 2015;Dasgupta 2021).
Freshwater ecosystems are one of the world's most vulnerable systems (Dudgeon et al. 2006;Khamis et al. 2014). While freshwater habitats occupy less than 1% of the planet, they contain 10% of all known species, 40% of fish diversity, and 30% of all vertebrate species (Dudgeon et al. 2006;Strayer and Dudgeon 2010). A relevant example is the case of wetlands, estuaries and coastal ecosystems, suffering from impacts that affect the status of at least 50% of these environments (Zedler and Kercher 2005;Barbier et al. 2011;Díaz et al. 2019).
A major threat to aquatic ecosystems is associated with the deterioration of water bodies due to human activities that cause habitat degradation (Vörösmarty et al. 2010;Arthington 2012;Collen et al. 2014;Davis et al. 2015). Underground water systems, which support several freshwater ecosystems, are facing substantial overexploitation, which creates unprecedented pressures on these resources but also on linked aquatic ecosystems. Groundwater resources are the main reserves of liquid freshwater on Earth. These resources support the 43 percent of irrigation agriculture worldwide and are the main source of drinking water for 50 percent of the global population (UNESCO 2015). However, increasing water demands during the last decades have caused a sharp intensification in groundwater withdrawals ending up in a widespread decline in their reserves and in the quality of these resources. Aquifers depletion is considered as a global phenomenon, with estimations alerting than the 20 percent of aquifers worldwide are already overexploited, and with withdrawals substantially exceeding recharge rates in arid and semi-arid regions (Famiglietti 2014; Thomas and Famiglietti 2019).
Despite the existence of a significant range of ecosystems with varying dynamics between them, all of them require an acceptable state of water bodies. Heavily groundwater pumping has resulted in falling water tables and underground bodies contamination (e.g. salinization, chemical pollution, nutrient pollution) with associated impacts in natural streamflows, hydrological systems, and ecosystems (Konikow and Kendy 2005). Several studies have already reported how groundwater depletion is an important threat to many freshwater ecosystems. The impacts include sharp decreases in wetlands' flooded areas and drying springs (Esteban and Albiac 2011;Cooper et al. 2015), reductions in the abundance and type of wetlands' vegetation (Stromberg et al. 1996;Froend and Sommer 2010), land subsidence (Dinar et al. 2021), penetration of invasive species (Danielopol et al. 2003), or river flows declines (Glazer and Likens 2012) among others.
The aim of this paper is to assess an optimal management of groundwater withdrawals while preserving aquatic ecosystems that present regime shifts. Our contribution lies in the consideration of switches in the state function of an ecosystem when a critical water level ('tipping point') is reached. Despite linear ecosystem behavior was considered a good approximation of the performance of ecosystems (Esteban and Dinar 2016), several studies point out that nonlinear specifications including shifts, tipping points, and even hysteresis processes, better represent the behavior of ecosystems (Scheffer et al. 2001;Scheffer and Carpenter 2003). For simplicity, in this manuscript we assume that groundwater depletion presents two alternative stable states after a certain tipping point is reached (e.g., decrease in the groundwater table level). Despite the fact that hysteresis processes and several tipping points could be a more realistic assumption (Scheffer et al. 2012;Folke et al. 2004), there are some impediments for the correct empirical implementation of that hypotheses, such as the lack of knowledge on the relationship between groundwater dependent ecosystems (GDEs), the groundwater regime, and the safe limits of groundwater reserves to maintain ecosystems (Eamus and Froend 2007). Furthermore, some studies have also reported the existence of little evidence for multiple states and tipping points (Schröder et al. 2005).
Our results highlight how optimal water extractions and aquifer water level paths are significantly different when a switch in the ecosystem status functions is taken into account. In contrast with previous contributions, our results suggest that optimal water extractions present discontinuous patterns. These outcomes contribute to a better understanding of the relationship between groundwater and ecosystems, and can be helpful for the implementation of effective groundwater policies to maintain economic activities while protecting aquatic ecosystems.
The remainder of the paper is organized as follows. Section 2 presents the hydro-economic representation of an aquifer including the dynamics of linked aquatic ecosystems. The resolution of the optimal control problem is stated in Sect. 3 and Sect. 4 describes the study area. Section 5 numerically illustrates the main theoretical findings and sensitivity analyses are performed in Sect. 6. Finally, Sect. 7 sets out the conclusions.

The Model
Depletion and quality deterioration of groundwater bodies has prompted a general interest in the analysis of an efficient management of these resources (see the reviews by Koundouri (2004) and Koundouri et al. (2017)). However, the effective management and control of groundwater resources is still far from being achieved.
Several studies have included ecosystems as relevant elements in groundwater management (e.g., Roumasset and Wada 2013;Gutrich et al. 2016;Esteban and Dinar 2016;Pongkijvorasin et al. 2018;Pereau et al. 2019). However, to the best of our knowledge none of them have specifically included the possibility of tipping points involving regime shifts in the ecological systems. de Frutos et al. (2019) analyze the occurrence of shocks and regime shifts in a groundwater resource, however, each shock is based on a sudden change in the dynamics of the resource itself. Our motivation in this paper is to include some notions of the ecological literature that has already reported how when a 'critical' habitat condition (tipping point or threshold) is reached, some systems switch from one regime to another, with important losses in the ecosystem state (Scheffer et al. 2001). Small modifications in the environmental conditions may cause dramatic shifts in the ecosystem ecological state (Chaparro-Pedraza and Roos 2020).
The economic literature is quite extensive in the analysis of managing natural resources in presence of regime shifts (e.g., Folke et al. 2004;Brozovic and Schlenker 2011;de Zeeuw and Zemel 2012;Crépin et al. 2012;di Maria et al. 2012;Werners et al. 2013;Ren and Polasky 2014;de Zeeuw 2014;Zemel 2014, 2017;Polasky et al. 2014;Lade et al. 2015), including the case of groundwater resources Zemel 1995, 2004;de Frutos et al. 2014 and. In general, most of the studies are mathematical and suggest 1 3 how optimal management patterns largely depend on the shifts in ecosystem. In this study we contribute with this literature by developing an optimal management model of an aquifer, linked with an aquatic ecosystem. This ecosystem presents a regime shift, once a water table threshold is reached, that originates a change in the structure and functioning of the ecosystem.

Hydro-economic Model
We take a typical single-cell aquifer with a flat bottom and a fixed natural recharge (R) that is used for human purposes (W t ) . The depth of the aquifer is represented by the difference between the surface level (S L ) and the aquifer water table level (H t ) . The groundwater dynamics (Ḣ) is a function of the natural recharge (R) minus net water extractions for human consumption (1 − ) ⋅ W t , with being a proportion of water that infiltrates the aquifer as a return flow from human extractions.
where AS represents the total available water (aquifer area multiplied by the storativity coefficient).
A social planner aims to maximize the net present value of future benefits streams from using the aquifer. Social benefits consist of both the private profits from human activities ( B W t , H t ) and the economic value of the goods and services from aquatic ecosystems ( E H t ). The functioning of aquatic ecosystems depends on the existence of an 'acceptable' groundwater level to guarantee their survival. So, the social planner problem model can be stated as follows: with W t being the total groundwater extractions (control variable) and H t being the groundwater table level (state variable). The social discount rate is denoted by r. This maximization is constrained by the dynamics of the resource (Eq. (1)), an the initial condition of the water table level H(0) = H 0 .
Private benefits from groundwater extractions are represented by a quadratic revenue function minus pumping costs (Koundouri et al. 2017). Traditionally, these withdrawals are related with extractions for irrigation activities that are the main groundwater user worldwide (Siebert et al. 2010).
Private net revenues from groundwater extractions ( Re W t ) are derived from a linear specification of a groundwater demand function W = g + k ⋅ P , where P is the groundwater price, and g and k are equation parameters ( g > 0 and k < 0 ). Total costs of extractions ( C W t , H t ) depend on both the amount of withdrawals and the depth of the groundwater level. C 0 and C 1 are respectively fixed costs and marginal costs from groundwater pumping (Gisser and Sanchez 1980).

Ecosystem Health Function
As reported in the literature, ecosystems' deterioration is non-linear and presents potential regime shifts when reaching certain habitat conditions (Scheffer et al. 2001(Scheffer et al. , 2012Beisner et al. 2003;Suding et al. 2004;Heffernan 2008;Utz et al. 2008;Crépin et al. 2012). Once a tipping point is reached, the ecosystem status switches to a different regime (Scheffer et al. 2001;Hughes et al. 2013). In this paper, we follow the representation by Scheffer and Carpenter (2003) on how the state of ecosystems changes under variations in the external conditions. Fig. 1 shows how ecosystems can respond smoothly to alternations in their habitats (a) or, by contrast, they can change abruptly when habitat conditions approach a critical level (b) with even having different stable states and several critical levels (c).
The ecosystem state in this study is represented by a function ( EC(S L − H t )) that links the health 1 of the ecosystem with the level of water into the aquifer (aquifer depletion). This function represents how decreases in the water table level (S L − H t ) affect the functioning of dependent aquatic ecosystems. When the aquifer is full (S L = H t ) the ecosystem maintains its pristine state; however, as the water table decreases, the health of the ecosystem declines (see Fig. 1b). In the case of groundwater dependent ecosystems, some studies state that aquifers depletion progressively impacts in the functioning and adaptation of ecosystems . However, 'if the stress is prolonged or extreme, these adaptations become inadequate and result in populations progressively declining, and a shift in the composition and function of ecosystems' (Rohde et al. 2017, p. 296). In this study we assume that the deterioration of habitat conditions (groundwater depletion) cause a profound change, non-linear relationship, in the ecosystem dynamic that presents two 'stable states' once a unique tipping point is exceeded.
The stated ecosystem function is described in Fig. 2 and shows how when the aquifer is full the ecosystem presents is maximum health ( 1 ) . The parameter 1 represents the ecosystem's best state (or health) when the aquifer has not been altered. On the other hand, the parameter 2 is the tipping point where the ecosystem switches to a different regime or state. Following the ecology literature, regime switches take place because external forces (e.g., habitat conditions) drive ecological conditions beyond a tipping point or critical threshold. In this model, the critical threshold (H C ) represents the water table level at which habitat conditions hamper the good status and viability of the ecosystem. Critical aquifer depletion depth (d c ) is defined as the difference between the surface level (S L ) and the critical water table level (H c ). Once the critical threshold is reached the ecosystem change to a different regime where their ecological status presents a notable deterioration. Following Scheffer et al. (2012), in the second state the marginal deterioration of the ecosystem is smaller because the remaining ecosystem presents a different ecological structure not as rich both in productivity and diversity as in the first state (Rohde et al. 2017).
The mathematical expression of the ecosystem health function EC S L − H t shown in Fig. 2 is as follows: 2 where 1 and 2 are positive parabola parameters, (S L − H t ) represents the aquifer depletion, H C is the critical threshold (critical water table level), and d c is the aquifer depletion depth at the critical threshold (S L − H c ) . We assume that at the critical threshold the ecosystem function has the same value both at right and left of the function. Finally, we assume that when the aquifer is totally depleted (H t = H bottom ) the ecosystem is extinguished. Even though several functional forms can characterize the representation of Fig. 1b, we state a parabola as a functional form of the ecosystem behavior to obtain mathematical solutions for the optimal control problem. In order to find theoretical solutions for the groundwater model presented in sub-Sect. 2.1, it was necessary to establish a quadratic function for the ecological functioning (Eq. (4)). Ecosystem State 2 We have defined the ecosystem health as function of the aquifer depletion (Fig. 2). The mathematical expression of Eq. (4) allows us to represent the ecosystem behavior based on biological literature (Fig. 1b). Equivalently, it could be possible to relate the health of the ecosystem with the water table level, and then an analogous equation for the ecosystem health function is: Finally, the economic value of the ecosystem depends on the amount of goods and services that the ecosystem provides. The economic value of the ecosystem is a linear function of the ecosystem state (EC): where is a constant parameter relating the ecosystem state and the aggregate economic value of the ecosystem services.

Two-stage Optimal Control Problem
Assuming the equations formulated in Sect. 2, the social planner's maximization becomes a 'two-stage' optimal control problem. When abrupt changes occur and the objective function shifts, the problem requires a resolution by phases ('multi-stage' dynamic optimization problems). 'An optimal multiprocess problem is a dynamic optimization problem involving a collection of control systems coupled through constraints in the endpoint of the state tra-jectories…' (Babad 1995, p. 530). We define this two-stage optimal control problem as the maximization of social welfare (SW), and we rewrite Eq. (2) as follows: with t c being the time at which the critical threshold is reached. The time t c is also one of the variables to be determined.
The 'two-stage' optimal control problem, corresponding with the two ecosystem's regimes, is solved by imposing a sequence of two phases with two Pontryagin problems (Tomiyama 1985;Tomiyama and Rossana 1989;Makris 2001;Boucekkine et al. 2004;Aisa et al. 2007). The 'two-stage' method proceeds by solving the two sub-problems backwards, with then the second optimal control sub-problem ( SP 2 ) resolved first: with W 2 and H 2 being the solutions of water extractions and water table level under the Pontryagin problem corresponding to sub-problem 2 ( SP 2 ). 3 The associated Hamiltonian is: where 2 is the costate variable of sub-problem 2. By solving the first order conditions we obtain the optimal results of this problem ( SP * 2 W * 2 , H * 2 , t C ), and the optimal values of the state and control variables ( W * 2 and H * 2 ) 4 : is the negative root of the polynomial equation from the system of two differential equations (see Appendix 1, Eq. (29)). Finally, by imposing the initial conditions in t c the coefficient B is obtained: 5 Having obtained the optimal value of sub-problem 2 ( SP 2 ), the second phase requires optimizing the Pontryagin problem while corresponding to the first sub-problem ( SP 1 ): 4 A complete resolution of sub-problem 2 (phase 1) is given in Appendix 1. 5 See the detailed resolution in the Appendix 1.
where W 1 and H 1 are the solutions under the sub-problem 1 ( SP 1 ).
The Hamiltonian associated with the Pontryagin problem corresponding to sub-problem 1 is: with 1 being the costate variable of sub-problem 1.
By solving the first order conditions, and imposing both the continuity condition * and the matching condition , we obtain the optimal paths for groundwater withdrawals (W * 1 ) and the water table (H * 1 ) under sub-problem 1. Additionally, we also obtain the optimal time (t * c ) for reaching H c . 6 The parameters CA and CB are constants to be determined with the initial conditions, Finally, y 1 and y 2 are the roots of the polynomial equation established to solve sub-problem 1. 7 The value of the constants CA and CB can be expressed as: The optimal paths for the water extraction (W * ) and the water table level (H * ) depend on the solutions of the two sub-problems ( SP 1 and SP 2 ) included in Eqs. (7) and (12). (13) Consistent with previous analysis (Ren and Polasky 2014;Crépin et al. 2012), our outcomes suggest that the inclusion of ecosystem regime shifts will condition groundwater management strategies when the objective is both, to maintain human activity and to preserve aquatic ecosystems. The theoretical outcomes show that both optimal variables (extractions and groundwater level) are conditioned by the ecosystem function parameters 1 and 2 (Eq. (4) and Fig. 2). Furthermore, the tipping point or critical threshold (H c ) is also a relevant element in the optimal equations.
The theoretical outcomes clearly indicate the relevance of the ecosystem dynamics (shifts and tipping point) in groundwater management. However, the complexity of the theoretical results, with several interactions between the model parameters, makes it difficult to assess the size and impact of these elements on the optimal outcomes. To illustrate the influence of these parameters on the results, we have performed a numerical analysis applied to the Eastern la Mancha aquifer (Spain).

Study Area
The Eastern La Mancha aquifer is located in Southeastern Spain within the upper Jucar river basin. This aquifer is the largest groundwater reservoir in the Iberian Peninsula with an extension of 7300 km 2 . The interest in the aquifer stems from the fact that it is a unique case of a large aquifer in arid and semiarid regions being managed towards sustainability, due to the success of the collective action engaged by stakeholders (Esteban and Albiac 2011). This is the reason for this aquifer being an international benchmark that could be used like an example for the analysis.
The expansion of irrigation in the Eastern La Mancha aquifer between the 1970s and 1990s caused a substantial decline in the aquifer's water table which led to a reduction in the streamflows to the Jucar basin, and produced important ecosystem damages. The aquifer water table has dropped about 80 m resulting in large storage depletion, fluctuating around 2500 Mm 3 . The aquifer is linked to the Jucar River stream, and it used to feed the Jucar River with about 150 Mm 3 /year in the 1980s but the feeding has fallen below 40 Mm 3 /year at present (Pérez-Martín et al. 2014). The consequence is that the lower Jucar is undergoing severe problems of low flows, water-quality degradation, and damages to the Albufera wetland, which is the most important aquatic ecosystem in the basin.
The Albufera wetland is a freshwater lagoon with an area covering 2430 ha and supporting very rich aquatic ecosystems. It was declared National Park in 1986, included among the international important wetlands by the Ramsar convention in 1990, and recognized as special area for bird protection (ZEPA) since 1991 (Soria 2006). The Albufera is located in the lower Jucar basin beside the Mediterranean Sea, and it is mainly fed by irrigation return flows from the Jucar River (Kahil et al. 2016). 'The case of the Albufera illustrates the complexity of water management in some Mediterranean wetlands where the connection between irrigation systems and palustrine ecosystems creates an amalgam of contrasted visions, goals, functions and values, generating enormous difficulties to define shared strategies for sustainability' (Jégou and Sanchis-Ibor 2019, p. 504). The reduction in the contribution of the inflows of the Eastern la Mancha aquifer into the Jucar River are dwindling the environmental flows in many parts of the basin, and especially in the lower Jucar River resulting in severe environmental damages to the Albufera ecosystems.
The institutional developments in the Eastern La Mancha aquifer leading to collective action started when farmers became aware of the problems from aquifer depletion and responded by creating a water-user association aimed to jointly manage the aquifer. This response was driven by the call for control of extractions by the basin authority with the strong support of downstream farmers, the threats by the basin authority of not issuing water rights to pumping farmers, and the increase of pumping costs because of the falling water table. These facts have resulted in a progressive reduction in extractions of 100 Mm 3 after 2000, down to the recharge level (Table 1).
In Table 2, we show the economic and hydrologic parameters of the Eastern la Mancha aquifer used for the numerical simulation Albiac 2011 and. Parameter values for both the groundwater demand and supply functions have been estimated from Esteban (2010). We have actualized these values with recent information regarding current average prices of irrigation in the region, the electricity costs (C 1 ) and the fixed costs of operation and maintenance of the wells (C 0 ) . Hydrological parameters are provided by Esteban and Albiac (2011).
For the ecosystem function and due to the lack of knowledge on specific groundwater ecosystem dynamics in the study area, we have approximated the parameters of this function ( 1 and 2 ). We assume that when the aquifer is full the ecosystem's heath is the best ( 1 = 1) and then, the economic value of this ecosystem is the highest. Furthermore, because the lack of knowledge on specific values for possible tipping points in ecosystems, we assume that the ecosystem' shift occurs when the health of the ecosystem has fallen by half ( 2 = 0.5) . To support our results, we have conducted a sensitivity analysis to evaluate the influence of these parameters on the optimal outcomes. Regarding the critical threshold (H c ) , we consider that it takes place after a groundwater depth of 40 m, which means to assume an additional depletion of 10 m compared with current situation (H 0 ) . Detailed information on the relationship between ecosystems health and water regimes, as groundwater depth, is still lacking and existing estimations are very different between regions, climatic variables, hydrology, geology, soil, and vegetation Huang et al. 2019). However, some studies report that groundwater depletion threatens vegetation species because plants in the riparian zone that grew because of close proximity of the water table may not survive as depth to water increases (e.g., Stromberg et al. 1996;Froend and Sommer 2010). Some estimation suggests 20 m of groundwater depth as a critical threshold for the survival of selected vegetation species (Chávez et al. 2016 Figure 3 illustrates the optimal trends of water pumping and the water table in the Eastern La Mancha when comparing two scenarios: (1) baseline scenario where ecosystems are ignored, and (2) scenario where ecosystem dynamics are taken into account (ecosystem scenario). The outcomes in Fig. 3 show the existence of differences under both scenarios, by considering and not considering the ecosystem dynamics, especially by the optimal pumping rates. Under the scenario where no ecosystems dynamics are considered (dotted line), groundwater withdrawals decrease until reaching the steady state. Optimal paths of groundwater withdrawals show the 'trade-off' between maintaining private benefits from groundwater pumping together with the preservation of the resource (extraction costs By contrast, when ecosystem dynamics are taken into account in the model, results show an important decrease in the initial level of extractions. Optimal withdrawals should be more than 20 Mm 3 lower compared with the baseline scenario (370 Mm 3 ). Furthermore, optimal paths of extractions show a declining trajectory for the first 60 years; and from then on, optimal extractions sharply increase until the ecosystem threshold (tipping point) is reached. Once the critical water table level is exceeded, optimal withdrawals are slightly higher up, compare with the baseline scenario, to approach the steady state.

Empirical Application: Results
Lower groundwater withdrawals, when environmental objectives are included in the model, is an expected result already reported in the literature (Pereau et al. 2019;Esteban and Dinar 2016). However, our results show a decreasing pattern followed by a sharp increase in extractions. The explanation of this unexpected pattern is related with the low contribution of the ecosystem to the social welfare. In this paper, we have approximated the groundwater-related ecosystems as some native vegetation that is directly linked with the depth of the aquifer water table. We assumed a conservative valuation of the ecosystem and a moderate economic impact of the ecosystem in the social welfare. So, after a 'short' period of time, private benefits from groundwater pumping largely exceed the benefits from the ecosystem and the model sacrifices ecosystems to benefit private extractions. With the selected ecosystem valuation and discount rate the preservation of the ecosystem is essential during almost 60 years. From then on, private profits from groundwater extractions are preferred and the ecosystem is pushed to their deterioration (second 'stable state'). An example of this, a priori surprising result, is the fact that the model described have been constructed using economic values of the ecosystem that are not "extreme". We have allowed the model to always reach the tipping point and causing a significant deterioration in ecosystems. However, when imposing large values of some relevant parameters, namely H c or the economic value of the ecosystem ( ) , the results show that the system would never reach the second sub-problem (SP2). An example of this result is represented in Fig. 4. This figure shows how large valuable ecosystems require an optimal water table level of 677 masl. This scenario is attained when assuming and economic values of the ecosystem goods and services of 40 million €. In this case, the optimization problem never reaches the second sub-problem ( SP 2 ) or the critical threshold. This outcome suggests that valuable ecosystems that provide several commodities to society need to be protected in order to preserve their ecological status. Under this scenario, government interventions are required with policies able to reduce groundwater withdrawals in order to avoid ecosystems deterioration.

Sensitivity Analyses and Model Extension
To corroborate and assess the impacts of different parameter values, we performed several sensitivity analyses with the coefficients related to the ecosystem function. The first sensitivity analysis conducted uses higher and lower values of the selected tipping point ( 2 ) . We have also simulated the impacts of imposing different levels of the critical threshold (H c ) . Aggregate outcomes of the main results of the sensitivity analyses are presented in Table 3. The results in Table 3 illustrate how the consideration of ecosystems dynamics will change the aquifer depletion levels and also social welfare.
Finally, we have included an additional simulation that incorporates the impact on the optimal variables when assuming that the shift in the ecosystem impacts on the aquifer's recharge rate. This scenario supposes an extension of the model and highlights the role of connecting the ecosystem with the hydrological behavior of the aquifer. Table 3 Main results from the scenarios analyzed * To estimate the aquifer depletion we establish a time horizon of 20 years ('short-run' result). This result allows a better understanding of the differences between scenarios and is more useful for policy makers

Scenario
Water  Figure 5 presents the scenario when higher and lower values of the tipping point ( 2 ) are stated. The tipping point is related with the ecosystem' resilience to habitat modifications due to represents the point at which the ecosystem health is being significantly altered. As expected, optimal patterns for groundwater extractions follow similar trends than those presented in Fig. 3. However, significant differences arise in both the time of reaching the critical threshold and the total volume of extractions (Table 3).
Results from Fig. 5 shows that when imposing a higher tipping point ( 2 = 3 5 = 0.6 ), the critical threshold is reached earlier in time (t * c = 93) . Furthermore, groundwater withdrawals are higher before reaching the critical threshold, but once exceeded are lower compared to the other scenarios. Finally, even though the peaks in extractions are similar in the three scenarios, when imposing a higher tipping point the peak is slightly lower. By contrast, with a lower tipping point ( 2 = 2 5 = 0.4) the critical threshold is reached after 258 years and the aquifer depletion is lower compared with the other two scenarios (Table 3). Under this scenario we also observed the highest peak in withdrawals, although the difference is small. These outcomes indicate that the resilience of an ecosystem, which is related with the greater ability to maintain their properties and processes, conditions the optimal trajectories for both optimal variables. Optimal trajectories indicate how highly resilient ecosystems (lower tipping point ( 2 = 0.4) ), should be protected for a prolonged period of time. It should be noted that this result is calculated by assuming the same economic value for high and low resilient ecosystems. The optimal model prioritizes an ecosystem that is able to contribute to the social welfare for a longer period of time. Figure 6 illustrates the optimal paths for both groundwater withdrawals and water level under different critical thresholds. In this case, we have selected upper (H c = 655) and lower (H c = 645) critical thresholds. 8 The results show that while optimal patterns follow similar trends, there are notable dissimilarities between aggregate water extractions and the time reaching the critical threshold t * c . A relevant outcome from this scenario is the significant difference in the groundwater withdrawals peaks depending on the tipping point. The upper the critical threshold, the greater the peak in extractions once this level is reached.

Sensitivity analysis with the critical threshold (H c )
A lower critical threshold means that the ecosystem displays a higher resilience to habitat alterations. The results of a sensitivity analysis imposing a lower critical threshold (H c = 645) show how the decrease in groundwater extractions holds for a prolonged period of time. As expected, the optimal time at which the critical threshold is reached (t * c ) occurs later. Finally, a relevant result is the fact that the peak in groundwater withdrawals, once the threshold is reached, is less severe compared with the other scenarios. By contrast, when imposing a higher critical threshold (H c = 655) , meaning that the ecosystem displays a lower resilience to habitat alterations, groundwater withdrawals decrease at a higher rate. Once the threshold is reached, at an earlier stage compared with the other scenarios, the peak in groundwater pumping is more significant. These results corroborate those obtained in the previous analysis and give weight to the argument that the optimal strategy is to protect higher resilient ecosystems with a greater capacity of providing goods and services to society. As noted above, these results depend on assuming the same economic value (similar provision of good and services) of low and high resilient ecosystems.

Model Extension: Change in the Recharge Rate
A final scenario simulates the case of assuming that the change in the ecosystem' function will affect the recharge capacity of the aquifer. Under this scenario, the groundwater dynamics stated in Eq. (1) becomes: where R 1 represents the recharge rate in the pristine situation (Table 2) and R 2 is the recharge rate once the critical threshold is exceeded. 'Vegetation is responsible for a large portion of the variation in recharge…' (Kim and Jackson 2011, pp. 9). Some studies suggest that the recharge rate increases when vegetation disappears due to lower evapotranspiration (Li et al 2018). So, to test this hypothesis we establish a higher recharge in the second steady state of the ecosystems. We simulate an increase of 2% and a 5%. 9 The new formulation of the groundwater hydrology modifies the theoretical results calculated in section 3. However, the modifications are minimal since the recharge is constant in the model and its modification does not alter the optimality conditions of the two-state optimal control problem. The effects on the model resolution when including Eq. 18 are presented in Appendix 2. Figure 7 shows the optimal trajectory of the water table and groundwater withdrawals when the shift in the ecosystem health function impacts in the recharge level of the aquifer. Under a higher recharge rate in the second state, based on the assumption that plant's populations are reduced and thereby the evapotranspiration will be lower, the critical threshold is reached later in time compare to the baseline. The results from this simulation demonstrate that changes in the recharge rate modify the time at which the tipping point is approached but also the total aquifer depletion. The inclusion of this effect intensifies the impacts of including the dynamic aspects in groundwater management. It is important to note that for simplicity in the model (Sect. 2) we assume that groundwater withdrawals for ecosystems (environmental flows) are marginal. A potential sophistication of the model could be the inclusion of ecological withdrawals as stated by Pereau et al. (2019).

Conclusions and Discussion
'Nature across most of the globe has now been significantly altered by multiple human drivers, with the great majority of indicators of ecosystems and biodiversity showing rapid decline' (Díaz et al. 2019, p. 4). Large decays in ecosystems during the last few decades, Fig. 7 Optimal paths of groundwater withdrawals and water table levels when the shifts in the ecosystem health alter recharge rates. Note: Solid line shows the baseline scenario. Dotted line represents the simulation with a 2% increase in the recharge rate and broken line represents the simulation with a 5% increase in the recharge rate 1 3 especially freshwater habitats, have increased awareness of the need to protect these systems. One of the main factors affecting freshwater ecosystems is groundwater depletion. The pressures placed on these hydrological systems with which ecosystems interact, together with the important alterations in ecosystem habitats and the growing exploitation of some ecosystem commodities, is threatening the ecological status of many ecosystems. Despite the fact that the study of groundwater resources has been fully addressed in the literature, the relevance of ecosystems and how they affect groundwater management is still pending.
In this paper, we present a groundwater model where aquatic ecosystems, which depend on the status of groundwater resources, are incorporated in the problem. Our contribution lies in the consideration of ecosystems' dynamics, meaning that there is a shift in the state function of the ecosystems once a critical groundwater level is reached. Traditionally, simple and smooth ecosystem behavior was considered as a good approximation of the performance of ecosystems to habitat changes (Esteban and Dinar 2016). However, several studies suggest that nonlinear specifications with shifts are better approximation to the behavior of ecosystems (Scheffer et al. 2001;Crépin et al. 2012).
Results from our theoretical model and the numerical illustration suggest that current groundwater management strategies are not optimal due to the fact that they do not incorporate ecosystem dynamics. The outcomes from our analysis highlight that optimal groundwater management is affected by ecosystem dynamics and the emergence of critical habitat thresholds or tipping points. As expected, when ecosystems are taken into account, aggregate optimal extractions must be lower and the optimal water table level must be higher compared to a situation without ecosystems. However, when including a shift in the ecosystem state function, optimal water extractions follow a non-continuous trend. Before reaching the critical threshold, extractions are lower compared with the scenario where ecosystems are not included. In order to maintain the ecosystem's habitat and a 'valuable' ecosystem, extractions need to be controlled. However, as the health of the ecosystem slowly deteriorates, and the water table reaches the critical threshold, water extractions suddenly increase to higher levels, even exceeding the levels under the scenario where ecosystems were not accounted for. Additionally, very valuable or low-resilient ecosystems could lead to the outcome of never reaching the critical threshold. To avoid the critical level, a policy to control groundwater extractions is needed in order to prevent large aquifer depletion levels that destroy linked-ecosystems. A relevant result in terms of policy is the fact that high resilient ecosystems should have a greater protection due to their larger capacity to provide goods and services to society. It is important to note that, this finding is obtained even though the model assumes that the economic value of ecosystems is the same for high and low resilient ecosystems.
The policy outcomes indicate that managers should carefully plan groundwater strategies to maintain human needs while protecting ecosystems. In order to achieve an efficient allocation of groundwater resources, there is a need for better knowledge on the GDEs and their functions and relationships with groundwater resources. Additionally, ecosystems will also affect some hydrological variables, such as aquifer recharge. The main political implication of our results is that policy regulations are mandatory if a social planner internalizes the protection of aquatic ecosystems. This is indeed an obvious outcome that has already been stated in the literature (Esteban and Albiac 2011;Roumasset and Wada 2013;Gutrich et al. 2016;Esteban and Dinar 2016;Pongkijvorasin et al. 2018;Pereau et al. 2019). However, our findings involve new implications based on the inclusion of ecosystem dynamics in groundwater management. Variables such as the critical threshold or the resilience of the ecosystem are relevant variables to take into account in groundwater management.
The empirical application deals with the Eastern La Mancha Aquifer, which is a unique case of a large aquifer in arid and semiarid regions being managed towards sustainability. This is a compelling reason for using this aquifer as an international benchmark, setting an example of successful collective action. The results obtained in the Eastern La Mancha Aquifer can be useful for policymakers when implementing regulations to control groundwater withdrawals. First of all, water policies require an accurate knowledge of the ecosystems linked with the water bodies and the dynamics of these systems. Second, better information on the interactions between groundwater bodies and ecosystem processes is needed. Third, managers need to know how different disturbances affect ecosystems, and also the resilience of the ecosystems to face these disruptions. And finally, the main challenge for institutions dealing with the environment and policymakers worldwide is to determine which is the 'acceptable' state of water bodies and how to entice the collective action of stakeholders.
The solution to this homogeneous differential equation is: where the parameters A and B are constants to be determined with the model initial conditions. x 1 , x 2 are respectively the roots of the polynomial equation (x 2 − r ⋅ x + n ⋅ m = 0).
By integrating Eq. 30 we obtain the optimal solution of the water table problem ( H 2 (t) ) under the first sub-problem (SP 2 ): To determine A , we use the transversality condition (Eq. (24)). This condition can only be satisfied when A = 0 (similar to the result by Gisser and Sanchez (1980)). 10 The coefficient B is determined by using the sub-problem initial condition H 2 t c = H c . The value of B can be expressed as: The optimal solutions of the optimal extractions and water table level under the social planner problem (SP2) are: Having the optimal values for both the state and the control variables, we proceed to evaluate the value of the functional defined in Eq. (19): It is demonstrated that it is a convergent integral because x 2 < 0 and r > 0 . Solving the integral and the limit we obtain the optimal value of the sub-problem, SP * 2 W * 2 t c , H * 2 t c , t c . The second step of the 'two-stage' maximization method involves the resolution of the second sub-problem (SP 1 ) with knowing the optimal solution of the first sub-problem (SP 2 ). The first sub-problem (SP 1 ) becomes: (30) W 2 − r ⋅Ẇ 2 + n ⋅ m ⋅ W 2 = 0 (31) W 2 (t) = A ⋅ e tx 1 + B ⋅ e tx 2 (32)