Marie-Jeanne Perrot-Minnot

Frank Cézilly

Parasite-Induced Behavioral Change: Mechanisms

Keywords: Crustaceans, Ethopharmacology, Insects, Neuromodulator, Parasite, Proteome, Rodents, Serotonin, Transcriptome

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Parasites can manipulate the behavior of their hosts to their own benefitthis is what evolutionary parasitology studies tell us. But let us go a step further and take up the challenge raised by these manipulative parasites messing with the brains of their hosts and giving our own brains a serious puzzle. How can a so-called 'simple' (not 'regressed') parasite hijack the behavior of its host, which in some instances might be a so-called 'higher' vertebrate? Is there anything like a 'manipulative molecule' secreted by the parasite to directly target its host'sC N Sa n ds p e c i fically modulate the behaviors affecting transmission success? Or does manipulation come as a fortuitous side-effect of the infection on the host immune system and metabolism?

Despite the growing number of studies reporting on behavioral manipulation by parasites, the proximate mechanisms underlying this phenomenon have been investigated in only a few of them (Table 1), and no mechanism has been completely elucidated. This article aims at reviewing these few cases, where the mechanisms of parasite manipulation have been investigated. However, we use here a broader approach, looking for a causal connection between altered host behavior and the modulation of gene expression in both the host and the parasite. The analysis is focused on the mechanisms underlying changes in behavior that increase parasite transmission success (strictly speaking, parasite-induced behavioral manipulation). Mechanisms associated with disease-related behavioral disorders, such as immune-generated alteration of the CNS, will not be addressed here.

The phenomenon of parasite manipulation can be fully understood only if the demonstration of a fitness gain for the parasite (ultimate cause) is coupled to the identification of the mechanisms underlying the observed behavioral changes (proximate causes). Understanding proximate causes of manipulation will contribute to our evolutionary analysis in two ways: (1) it will help evaluate the costs a parasite pays to invest in manipulation and whether these costs are shared with investment in parasite survival (i.e., defense against the immune system) and (2) it will reveal how complex and specific the manipulation process is. These two criteria are currently acknowledged as important in assessing the adaptive significance of manipulation.

The mechanisms involved in parasite manipulation have been explored since the pioneering work on rodents infected with Toxoplasma gondii and on the amphipod Gammarus lacustris infected with the acanthocephalan bird parasite Polymorphus paradoxus. Since then, several studies have attempted to identify the changes in host neurophysiology or gene expression associated with parasite manipulation (Table 1; [START_REF] Tain | Altered host behaviour and brain serotonergic activity caused by acanthocephalans: Evidence for specificity[END_REF]. The expected complexity of the interactive network connecting a host and its manipulative parasite comes from the modulatory connections between the neuronal, hormonal, and immune systems of the host. The investigation of proximate mechanisms therefore relies on an integrative approach combining behavioral ecology, neurophysiology, pharmacology, molecular biology, and biochemistry. (5-HT) [START_REF] Helluy | Serotonin, octopamine, and the clinging behavior induced by the parasite Polymorphus paradoxus (Acanthocephala) in Gammarus lacustris (Crustacea)[END_REF] Positive phototactism [START_REF] Maynard | Gammarus lacustris harboring Polymorphus paradoxus show altered patterns of serotonin-like immunoreactivity[END_REF] Source: a: NO is a neuromodulator involved in memory, neuronal development as well as immune defense.

b: ATP can also play a role as neuromodulator. c: Inducing apoptosis in the brain can alter the chemical signals in the brain (inflammatory immune responses), as reported in vertebrates infected with certain viruses or protozoans.

From Phenotypic Behavioral Changes to Altered Gene Expression

Changes in host behavior following infection are not necessarily profitable to the parasite. They may actually benefit the host through compensating for the effect of infection or getting rid of the parasite. Alternatively, they can be pathological side-effects, with no benefits for the host or for the parasite. When beneficial to the parasite, changes in host behavior can result from a combination of direct and indirect effects of a parasite on its host's CNS. The most likely indirect effect relies on the connection between the neuronal and immune systems: the host's immunological response to infection can be involved in changing the host's behavior into a behavior that favors parasitic transmission. Therefore, the methods used to investigate mechanisms of parasite-induced behavioral changes must not only identify the biochemical or physiological changes in manipulated hosts, but also demonstrate that these changes are indeed the proximate cause of behavioral manipulation.

More precisely, we have to identify the following

(1) The functional connection between an altered behavior and the corresponding genes expressed in the host.

(2) The parasite's biochemical signal (in the excretory/secretory (E/S) parasite products) targeting host's genes, whether it corresponds to 'manipulative molecules,' or molecules with a broader spectrum, including physiological targets.

(3) The causal link between some of those genes and the direct target of the parasite's E/S products.

To that end, several complementary approaches are possible

• The exploration of specific neurophysiological pathways by means of ethopharmacology and techniques used on candidate proteins (immunohistochemistry or immunocytochemistry, HPLC-ED, etc.). This approach can establish a functional link between a neuromodulatory pathway and the observed altered behavior, without establishing how the parasite directly hijacks the neurophysiology of its host.

• The differential screening of the host proteome or transcriptome between infected-manipulated individuals and non- manipulated ones (uninfected and infected), to reveal proteins or mRNA associated with a manipulated phenotype (as the cause or the consequence of altered behavior and physiology).

• The proteomic analysis of parasite's E/S products followed by the identification of the biological fractions modulating host's behavior. Proteomic tools applied to the analysis of E/S products screen for molecules released by a manipulative parasite that could trigger the observed phenotypic changes.

The first two approaches must compare infected-manipulated individuals with nonmanipulated ones (infected by a nonmanipulative stage of parasite and uninfected), to specifically identify neurophysiological or biochemical changes associated with manipulation. Still, will these pathways or molecules in the host's repertoire be the direct target of parasite? The third approach is thus necessary to identify the initial molecular dialog setting up behavioral manipulation.

Several inferences emerge from the astonishing fact that parasites increase their own transmission success by taking control of their host's behavior. (1) The molecular cross-talk between a host and a parasite that results in fine-tuned phenotypic alterations is probably complex and intimate. (2) A parasite manipulating the behavioral flexibility of its host so that it performs the appropriate behavior likely uses either molecular mimicry (biochemical evolutionary convergence) or highly conserved molecules (phylogenetic inertia). The parasite may thereby control some of its host modulatory pathways by usurping signaling processes. (3) Changes in host behavior are often a mix of direct and indirect effects, and it may prove difficult to differentiate between the two. Investigations of proximate mechanisms involved in parasite manipulation must keep these inferences in mind.

Multidimensionality and Mechanisms of Parasite Manipulation

The capacity of a parasite to manipulate several behavioral and physiological traits together has been largely ignored in most empirical studies so far, although a review of studies on the same host-parasite systems shows that manipulative parasites generally modify more than a single dimension in the phenotype of their hosts. For instance, the acanthocephalan fish parasite Pomphorhynchus laevis reverses the photophobic behavior of its host Gammarus pulex and its antipredatory behavior in reaction to olfactory cues, and increases its activity and its drifting behavior. Several physiological changes have been reported as well in G. pulex infected with P. laevis, such as increased hemolymph protein titers (in particular haemocyanin), reduced O 2 consumption, increased glycogen content, fecundity reduction, and immunosuppression. In wild rats infected with the protozoan T. gondii, changes in activity and in motivational level in various contexts have been reported. T. gondii-infected rats were found to be significantly less neophobic toward food-related novel stimuli. In outdoor captive environment, they were more likely to be trapped than their uninfected conspecifics, and their propensity to approach a mildly fear-inducing object was higher than that of uninfected rats (reduced neophobia). Alteration of innate behavior (such as neophobia) extends to the reversal of antipredatory behavior from a strong aversion to a preference for cat-treated areas in infected rats. This 'fatal attraction' is expected to increase the chances of transmission of T. gondii to its feline definitive host. Such multidimensionality of manipulation makes sense from an ecological and evolutionary point of view: having the 'vehicle' host reaching the right place at the right time (through being predated by, or stinging, or biting the next host species in the cycle) probably involves several behaviors related to environmental sensing and microhabitat choice. In parasites with a direct life cycle, transmission by contact or wounding can be increased by modulating a number of social behaviors, such as aggression and exploration. Several cue-oriented behaviors are generally altered in infected invertebrates (among phototaxis, chemotaxis, rheotaxis or wind-evoked behavior, geotaxis, etc.) that together contribute to increased transmission success of the manipulative parasite.

Are these multiple dimensions of a manipulated phenotype functionally independent? Or do the proximate mechanisms of manipulation have 'pleitropic effects'? The best argument supporting the hypothesis of 'pleiotropic effects' lies in the functional connection between host's neuronal, immunological, and endocrinal/metabolic systems, be the host an invertebrate or a vertebrate. Because the very first conditions for a parasite to develop are to successfully establish in a host and exploit its energy reserves, some mechanisms must exist that allow the parasite to interact with its host's physiology, especially the host's immunity. As pointed out by several authors, the evolutionary transition leading to parasite manipulation may simply consist in an extension of the effect of the parasite on the immune system of its host to its neuronal system. Targeting diverse and flexible neuromodulatory pathways to induce adaptive behavioral change in its host would thereby be a small evolutionary step. The understanding of proximate mechanisms of parasite manipulation allows us to test this evolutionary and functional scenario.

Investigating Host's Neuromodulatory Pathways

Biogenic amines (serotonin, dopamine, octopamine among others) and other chemical signals such as neuropeptides or the gas nitric oxide (NO) play a neuromodulatory role in numerous sensory, motor, and endocrine functions, in both invertebrates and vertebrates. They modulate the behavioral or physiological responses of an organism to external information, according to its internal status. By 'manipulating' these neuromodulatory pathways in its host's CNS, a parasite could adjust its host's behavioral response to reflect the parasite's own interest. The E/S products of the parasite would thus be akin to the venom of several predators or parasitoid wasps manipulating the monoaminergic system of their hosts to improve prey handling and use.

Several studies have shown a major role of biogenic amines and neuropeptides in the physiological and behavioral alterations induced by parasites (Table 2). The 'candidate neuromodulatory pathway' approach to parasite manipulation targets simple tropisms or cue-oriented behaviors such as phototaxis, geotaxis, chimiotaxis, thermal gradient sensitivity (in biting or sucking vectors of warm blood animals), and reflectance (Table 2). In vertebrates, several viruses and protists increase their hosts' exploratory behavior or aggression, two behaviors suspected to enhance parasite transmission either by predation or by conspecific wounding/contact respectively. These behavioral effects have been related to changes in concentrations or receptor binding of amines (dopamine, serotonin) or opioids.

However, few studies have combined ethopharmacological analysis to biochemical techniques (immunohisto-or cytochemistry, western blot, and ELISA) to establish or invalidate the involvement of a neuromodulator in parasite manipulation of behavior. One pioneer ethopharmacological study investigated the role of several neuromodulators in the behavioral alterations induced by P. paradoxus in its intermediate host G. lacustris. Uninfected individuals injected with serotonin responded to mechanical stimulation by skimming to the water surface until clinging to floating material and exhibited positive photaxis, two behavioral mimics of amphipods infected with this parasite of dabbling ducks. Immunocytochemistry on the nerve cord of amphipods infected with P. paradoxus revealed an increase in the number of varicosities exhibiting serotonin-like immunoreactivity in the third thoracic ganglion. Serotonin was altered either in the amount or in the number of local storage and release sites along neural fibers in P. Paradoxus-infected amphipods, but not in G. lacustris infected with Polymorphus marilis, a parasite of diving ducks inducing positive phototactism but no escape response.

Exogenously supplied serotonin can mimic the effect of parasitism in other amphipod-acanthocephalan systems: G. pulex-P. laevis and Pomphorhynchus tereticollis. Injection of serotonin in uninfected G. pulex reversed their reaction to light, mimicking the positive phototactism of gammarids infected with these two fish parasites. The serotonergic activity in the brains of infectedmanipulated gammarids was significantly increased, compared to that of four controls: uninfected G. pulex, G. pulex infected with P. tereticollis but not manipulated, G. pulex infected with the bird acanthocephalan Polymorphus minutus (not altering phototaxis), and a nonmanipulated sympatric amphipod species Gammarus roeseli, infected with P. laevis (Fig. 1).

In Gammarus insensibilis infected by the cerebral trematode Microphallus papillorobustus, immunocytochemistry on brain has revealed the degeneration of discrete sets of serotonergic neurons: immunoreactivity to serotonin (5-hydroxytryptamine or 5-HT) was decreased in the optic neuropils but increased in the olfactory lobes. This imbalance in brain serotonergic activity is suspected to contribute to the behavioral alterations reported in this brackishwater amphipod species, in particular, positive geotactism and attraction to light. In vertebrates, several viral and protozoan parasites infecting the CNS of their rodent hosts alter neurochemical pathways in the brain. In the brains of infected mice and rats for instance, rabies virus decreases 5-HT and GABA neurotransmission, and T. gondii increases the concentration of dopamine and decreases the concentration of norepinephrine. These changes in neuromodulatory pathways may be linked to elevated aggression exhibited by infected rodents (and exploratory and fearless behavior in the case of T. gondii infected rats). These behavioral alterations presumably enhance the transmission of rabies virus by increased conspecific biting and of T. gondii by increased predation rate.

Although the exploration of these neurophysiological changes can provide evidence that a neuromodulator plays a key role in one or few behavioral dimensions of parasite manipulation, it also has several limitations. First, the neuromodulatory and signaling network is complex: several neuropeptides or amines may act together to modulate a given behavior, while a single neuromodulator may regulate several behaviors. If this may fit well with the multiple dimensions of parasite-induced alteration on host's phenotype, it makes the full understanding of the underlying neurophysiological process difficult. Second, showing a change in brain CNS does not establish a causal connection with the manipulative process.

Table 2

Review of studies suggesting the involvement of certain neuropeptides or biogenic amines in parasite-induced alteration of invertebrate hosts' behavior (used as intermediate hosts by trophically transmitted parasites or as food store and shelter by parasitoid larvae)

Host species

Parasite species -Injection of octopamine or blood from postemergence parasitized larvae -Injection mimics the decreased peristaltic activity in the foregut (related to decreased feeding) [START_REF] Miles | Octopamine mimics the effects of parasitism on the foregut of the tobacco hornworm Manduca sexta[END_REF] Periplaneta americana

Technics

Results

Screening the Host's Proteome and Transcriptome

In the few host-parasite systems to which it has been applied, the proteomic approach appeared sensitive enough to detect proteome differences between infected and noninfected hosts that can be attributed to the manipulative syndrome. Indeed, the use of comparative screening of whole proteome or transcriptome between infected hosts and uninfected hosts appears a powerful means to cope with the predicted complexity of proximate mechanisms involved in parasite manipulation, if several conditions are met (e.g., the quality of controls run, the access to database allowing protein identification, and other limitations listed here). Proteins or transcripts differentially produced and specifically associated with the manipulative process can be identified, if one compares manipulated hosts with uninfected and infected nonmanipulated hosts. The analysis of infected nonmanipulated hosts (i.e., usually containing a developmental stage of the parasite not infective to the next host) is an important control to run, to distinguish the proteins or transcripts specifically associated with the manipulative process from the ones produced in response to infection. Similarly, noninfected hosts exposed to the same environmental conditions as infected manipulated ones should be analyzed (in addition to noninfected ones in their natural environment) to distinguish the proteins or transcripts specifically associated with the manipulative process from the ones produced in response to the environmental changes associated with manipulation (for instance, living at the surface instead of the bottom of a body of water). The differences in brain proteome between infected-manipulated hosts and controls are either in the presence/absence, the quantity, or the posttraductional processing of certain proteins.

From the studies reviewed (Table 3) it seems that the alteration of the CNS is a common feature in the proteome of infected manipulated animals. In addition, key metabolic pathways are often perturbed, as well as proteins involved in cellular stress (HSP, other chaperones), immunomodulation, or oxidative damage. Alteration in energy metabolism in the brains of infected blood-feeding insects can be interpreted as a parasite strategy to manipulate vector-feeding behavior by inducing a nutritional stress. Interestingly, several proteins putatively involved in similar behavioral modifications in different host-parasite systems belong to the same family. For instance, differential expression of proteins from the family (CRAL-TRIO) implicated in the vision process has been found in the brain of the wood cricket N. sylvestris infected by the manipulative nematomorph P. tricuspidatus and in the brain of G. insensibilis infected with the bird trematode M. papillorobustus (Table 3). Such pattern is suggestive of a limited 'jeu des possibles,' with molecular convergence or conserved proximate mechanisms being the only way to alter the behavior of phylogenetically distant hosts such as an amphipod and an insect.

However, despite their power in investigating the molecular basis of parasite manipulation, proteomic studies have several limitations

• Some peptides might be undetected because of their low concentration, specific pI, or small size (for instance, neuropeptides).

• Most proteomic studies have focused on differential expression of proteins; however posttranslational modifications of proteins might be involved as well in the modulation of the host's phenotype induced by parasites, and their importance in interpretating transcriptional data has been emphasized by some authors.

• Protein identification becomes problematical in nonmodel organisms: it relies on cross-species protein identification, and therefore, on either highly conserved proteins or proteins known from organisms closely related to the host. For instance, the power of proteomic studies on Gammarus amphipods in revealing proteins specifically linked to the manipulative process was limited by the impossibility of identifying 27 out of 72 proteins spots differentially present or absent from the brain of G. insensibilis manipulated by M. papillorobustus and 60 out of 68 proteins spots differentially present or absent from the brain of G. pulex manipulated by P. minutus.

Microarrays provide an alternative to proteomics but have been used so far in a limited number of host-parasite systems (rodents, salmons, mosquitoes, and bees infected with Apicomplexa protists, helminths, or mites). Most studies have performed transcriptional profiling of infected and uninfected individuals or resistant and susceptible strains, to identify the immune and metabolic response to infection. In parasites, microarray profiles coupled to other techniques in functional genomics have been used to discover new drug targets, or to understand the genetic basis of drug resistance. To our knowledge, the method has not yet been used to screen for 'manipulative molecules,' possibly because of the lack of genomic data on the historically best model systems of behavioral manipulation (acanthocephalans or trematodes and their arthropod intermediate hosts). This could, however, be done on mosquitoes/Plasmodium (with altered behavior driven by olfactory cues), tsetse flies and Trypanosoma, and bee/Varroa systems. Indeed, in the later model, a recent microarray-based analysis reveals that resistant and susceptible strains of bee to Varroa differ more on olfaction-related genes than on immune genes.

The Search for the 'Molecules of Manipulation' in E/S Products of Manipulative Parasites Differential screening of infected-manipulated hosts and noninfected or infected and nonmanipulated hosts will not reveal which molecules are actually released by the parasite to induce changes in the proteome or transcriptome. Even if proteomic studies can reveal a protein specifically produced by a parasite during the manipulation process, such as the mimetic protein Wnt by hairworms infecting crickets, the demonstration of its role as a 'manipulative molecule' requires that its excretion in the host's hemolymph is established. Therefore, proteomic analysis of E/S products from manipulative parasites is necessary to identify the parasitic molecules initiating the manipulation process. The separation of biological fractions in the 'secretome' allows their biochemical identification and must be followed by the biological testing in vivo of their functional role (in the manipulative process). Such analysis is currently limited to parasites that can be maintained or cultivated in vitro. The proteomic identification of E/S parasite products also relies on ongoing sequencing (genomic and EST) projects. So far, it has been possible in several nematodes and in the trematode Schistosoma mansoni, in the prospect of designing vaccines or drugs. Given these limitations, no such analysis has yet been done to identify the E/S products involved in behavioral manipulation. Studies analyzing E/S products from worms have revealed the release of proteins involved in a diversity of functions, such as stress response proteins/chaperones, antioxidant enzymes, energy metabolism and structural/cytoskeletal proteins, immune evasion, protease inhibitors, and lipid binding. No doubt the proteomic analysis of 'secretome' is a promising route to the discovery of the 'manipulative molecules' (either directly or indirectly targeting the host's CNS).

Conclusion

The most tricky aspect in our quest for the mechanistic basis of parasite-induced behavioral alterations is a causality problem akin, at first sight, to 'the chicken or the egg' dilemma. Whenever one spots a change associated with the manipulative process/the manipulated phenotype (be a neurophysiological change, a gene or a protein differentially expressed or produced), is it the cause or a sideeffect of parasite manipulation? The reviewed studies addressing the mechanisms underlying parasite manipulation reveal that the alteration of host's CNS is a common feature of manipulated hosts. Given the diversity of host-parasite systems in which behavioral alterations have been reported, further studies should help in answering several key questions: (1) Have the same constraints on parasitic transmission led to a similar solution (evolutionary convergence of behavioral manipulation)? (2) Are similar behavioral alterations induced by the same biochemical tools (molecular convergence)? For instance, geotaxis is altered in several host-parasite systems. Ants parasitized with parasitic fungus of the genus Cordyceps and ants infected with the liver fluke Dicrocelium, as well as the caterpillar Mamestra brassicae infected with M. brassicae nuclear polyhedrosis virus, climb to the top of a plant. Several amphipods infected with various acanthocephalan species (G. pulex/P. minutes, G. lacustris/P. paradoxus, and G. insensibilis/M. papillorobustus) swim at the water surface. Chemosensing and possibly learning are also altered in several host-parasite systems. Changes in olfactory perception and/or learning have been reported in various host-parasite systems as different as rats infected with T. gondii, amphipods infected with acanthocephalans, or dipteran vectors carrying protozoans. Is molecular convergence at the heart of these similar changes? The understanding of the mechanisms underlying parasite manipulation will likely reveal how parasites play tricks on their hosts, guiding them round more or less directly, by using molecular mimics or conserved molecules that affect the host's flexible neuromodulatory network.

  jumping into it Brain proteome (host) and parasite proteomeBiron et al.

-

  GC-MS and HPLC-ED on venom -Pharmacological depletion of monoamines mimics venom-induced nonparalytic hypokinesia and reduced-escape response -Immunocytochemistry (anti-dopamine) -Dopamine present in the venom -Electrophysiology -Decrease in octopamine neurons activity in the thorax, modulated by input from descending neurons from the brain, themselves modulated by venom injection Source: a: NO is a neuromodulator involved in memory, neuronal development as well as immune defense.b: ATP can also play a role as neuromodulator. c: Inducing apoptosis in the brain can alter the chemical signals in the brain (inflammatory immune responses), as reported in vertebrates infected with certain viruses or protozoans.

Fig. 1

 1 Fig. 1 5-HT immunoreactivity (yellow) within the brains of uninfected Gammarus pulex, P. tereticollis-infected G. pulex, and P. minutus-infected G. pulex. Arrows show position of tritocerebrum giant neuron (TGN) cell body. No differences in brain anatomy from infected and uninfected individuals were observed. Bar shows 100 mm.

  See also: Cognition: Non-Elemental Learning in Invertebrates. Hormones and Behavior: Immune Systems and Sickness Behavior. Host-Parasite Interactions: Disease: Parasite-Modified Vector Behavior; Evolution of Parasite-Induced Behavioral Alterations; Intermediate Host Behavior; Propagule Behavior and Parasite Transmission. Methodology: Neuroethology: Methods. Neurons and Senses: Nematode Learning and Memory: Neuroethology; Parasitoid wasps: Neuroethology. Predator-Prey Interactions: Parasitoids.

Table 1

 1 Review of studies attempting to identify the changes in host's neurophysiology or proteome associated with parasite manipulation in invertebrates

	Method of investigation References	Brain proteome (host) and parasite proteome Biron et al.
	Manipulated behavior	Seeking water and jumping into it
	Parasite	Hairworm, Spinochordodes tellinii	(Nematomorph)
	Host	Grasshopper, Meconema thalassinum	(Orthopteran insect)


Amine Chekireb
email: chekireb@cerege.fr

Julio Goncalves
email: goncalves@cerege.fr

Hubert Stahn
email: hubert.stahn@univ-amu.fr

Agnes Tomini
email: agnes.tomini@univ-amu.fr

Private exploitation of the North-Western Sahara Aquifer System

Keywords: Hydro-economic model, private pumping, multi aquifer system, groundwater-dependant ecosystems, semi-arid region, simulation. JEL classification: C61, C62, Q15, Q25

We formulate a hydro-economic model of the North-Western Sahara Aquifer System (NWSAS) to assess the effects of intensive pumping on the groundwater stock and examine the subsequent consequences of aquifer depletion. This large system comprises multi-layer reservoirs with vertical exchanges, all exploited under open access properties. We first develop a theoretical model to account for relevant features of the NWSAS by introducing, in the standard Gisser-Sanchez model, a non-stationary demand and quadratic stock-dependent cost functions. In the second step, we calibrate parameters values using data from the NWSAS over 1955-2000. We finally simulate the time evolution of the aquifer system with exploitation under an open-access regime. We specifically examine time trajectories of the piezometric levels in the two reservoirs, the natural outlets, and the modification of water balances. We find that natural outlets of the two reservoirs might be totally dried before 2050.

Introduction

Groundwater use has increased significantly worldwide, especially in many arid and semiarid areas where surface water is scarce. Intensive exploitation, combined with poor management, has caused a substantial drop in groundwater levels, water quality degradation, and significant environmental impacts. Water scarcity affects local activities that depend on water extraction and users located within natural discharge areas. Increasing groundwater extraction can thus threaten water flow and dry up natural springs on which people and ecosystems are dependent. This study investigates the hydrological impacts of unregulated groundwater extraction in an aquifer system. We mainly concentrate on the case of the North-Western Sahara Aquifer System (NWSAS).

Groundwater exploitation has been extensively analyzed in resource economics. Since the seminal work of Gisser and Sánchez [START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF], many studies have, however, adopted a more interdisciplinary approach to account for hydrological characteristics of the aquifer, and interactions within the hydrological system, as well as feedbacks between the hydrological system and the economic activity (Bekchanov et al. [START_REF] Bekchanov | Systematic review of Water-Economy modeling applications[END_REF], Brouwer [START_REF] Brouwer | Integrated hydro-economic modelling: approaches, key issues and future research directions[END_REF], Harou et al. [26]). Hydroeconomic models (HEMs) have been used to assess and implement water policies at different scales. 1 These models allow to characterize specific features of different water systems -river basins (Cai et al. [START_REF] Cai | Physical and economic efficiency of water used in the river basin: implications for efficient water management[END_REF], Pulido-Velasquez et al. [START_REF] Pulido-Velasquez | Hydro-economic river basin modelling: Ther application of a holistic surfacegroundwater model to assess opportunity costs of water use in Spain[END_REF], or Torres et al. [START_REF] Torres | Economic impacts of regional water scarcity in the Sao Francisco River basin, Brazil: An application of a linked hydro-economic model[END_REF]), groundwater aquifers (Medellín-Azuara et al. [START_REF] Medellín-Azuara | Hydro-economic analysis of groundwater pumping for irrigated agriculture in California's Central Valley, USA[END_REF]) or surface-groundwater interactions (Kahil et al. [START_REF] Kahil | Hydroeconomic modeling with aquifer-river interactions to guide sustainable basin management[END_REF], Stahn and Tomini [START_REF] Stahn | Rainwater harvesting under endogenous capacity of storage: a solution to aquifer preservation?[END_REF] [START_REF] Stahn | On the environmental efficiency of water storage: The case of a conjunctive use of ground and rainwater[END_REF]) -and to investigate explicit issues, including droughts and water scarcity (Maneta et al. [START_REF] Maneta | A spatially distributed hydro-economic model to assess the effects of droughts ion land use, farm profits, and agricultural employement[END_REF]), water allocation policies (Booker and Young [START_REF] Booker | Modeling intrastate and interstate markets for Colorado river water resources[END_REF], Chatterjee et al. [START_REF] Chatterjee | The optimal provision of water for irrigation and hydropower[END_REF], Esteve et al. [START_REF] Esteve | A hydroeconomic model for the assessment of climate change impacts and adaptation in irrigated agriculture[END_REF]), nitrate pollution (Graveline and Rinaudo [START_REF] Graveline | Constructing scenarios of agrcultural diffuse pollution using an integrated hydro-economic modelling approach[END_REF], or Peña-Haro et al. [START_REF] Peña-Haro | Stochastic hydroeconomic modeling for optimal management of agricultural groundwater nitrate pollution under hydraulic conductivity uncertainty[END_REF]), or recently climate change (Connor et al. [START_REF] Connor | Impacts of climate change on lower Murray irrigation[END_REF], Esteve et al. [START_REF] Esteve | A hydroeconomic model for the assessment of climate change impacts and adaptation in irrigated agriculture[END_REF]). Among all these models, the Gisser-Sánchez's HEM continues to be widely used to analyze potential benefits from optimal management. Indeed, a large strand of the literature addresses the Gisser-Sánchez effect, stating that policy regulation entails limited benefits. This is evident in a number of studies that has reproduced the Gisser-Sánchez HEM (e.g., Allen and Gisser [START_REF] Allen | Competition versus Optimal Control in Groundwater Pumping when Demand is Nonlinear[END_REF], Brill and Burness [START_REF] Brill | Planning versus Competitive Rates of Groundwater Pumping[END_REF], Feirnermann and Knapp [START_REF] Feinerman | Benefits from Groundwater Management: Magnitude, Sensitivity, and Distribution[END_REF], Koundouri [START_REF] Koundouri | Potential for Groundwater Management: Gisser-Sánchez Effect Reconsidered[END_REF], Negri [START_REF] Negri | The Common Property Aquifer as a Differential Game[END_REF], Nieswiadomy [START_REF] Nieswiadomy | The Demand for Irrigation Water in the High Plains of Texas, 1957-1980[END_REF], Provencher and Burt [START_REF] Provencher | The externalities associated with the common property exploitation of groundwater[END_REF], Rubio and Casino [START_REF] Rubio | Competitive versus Efficient Extraction of a Common Property Resource: The Groundwater Case[END_REF], Tomini [START_REF] Tomini | Is the Gisser and Sánchez model too simple to discuss the economic relevance of groundwater management?[END_REF]). Such optimization-based HEM has been then used to tackle other issues as uncertainty (Knapp and Olson [START_REF] Knapp | The Economics of Conjunctive Groundwater Management with Stochastic Surface Supplies[END_REF], Gemma and Tsur [START_REF] Gemma | The Stabilization Value of Groundwater and Conjunctive Water Management under Uncertainty[END_REF], Tsur and Graham-Tomasi [START_REF] Tsur | The Buffer Value of Groundwater with Stochastic Surface Supplies[END_REF]), quality (Roseta-Palma [START_REF] Roseta-Palma | Groundwater management when Water Quality is Endogenous[END_REF], Hellegers et al. [START_REF] Hellegers | Dynamics of Agricultural Groundwater Extraction[END_REF]), multiple water sources (Burt [11], Provencher [START_REF] Provencher | Issues in the Conjunctive Use of Surface Water and Groundwater[END_REF], Tsur [START_REF] Tsur | The Economics of Conjunctive Ground and Surface Water Irrigation System : Basic Principles and Empirical Evidence from Southern California, in Decentralization and Coordination of Water Resource Management[END_REF], Stahn and Tomini [START_REF] Stahn | Rainwater harvesting under endogenous capacity of storage: a solution to aquifer preservation?[END_REF]- [START_REF] Stahn | On the conjunctive management of groundwater and rainwater[END_REF], and Roumasset and Wada [START_REF] Roumasset | Ordering Renewable Resources: Groundwater, Recycling, and Desalination, The B.E[END_REF]). A separate strand of the literature deals with water flows and emphasizes the importance of groundwater flows (e.g., Brozovic et al. [START_REF] Brozović | On the spatial nature of the groundwater pumping externality[END_REF], Pfeiffer and Lin [START_REF] Pfeiffer | Groundwater pumping and spatial externalities in agriculture[END_REF], Esteban and Albiac [START_REF] Esteban | The problem of sustainable groundwater management: the case of La Mancha aquifers, Spain[END_REF], Guilfoos et al. [START_REF] Guilfoos | Groundwater management: The effect of water flows on welfare gains[END_REF], Perreau et al. [START_REF] Perreau | Optimality versus viability in groundwater management with environmental flows[END_REF] or Reinelt [START_REF] Reinelt | Spatial-dynamic seawater intrusion and pumping cost externalities in a confined aquifer[END_REF]). Guilfoos et al. [START_REF] Guilfoos | Groundwater management: The effect of water flows on welfare gains[END_REF] explicitly analyze the role of public policy of a multi-cell aquifer and outline the spatial effects on welfare gains with an application to Kern County (California). Esteban and Albiac [START_REF] Esteban | The problem of sustainable groundwater management: the case of La Mancha aquifers, Spain[END_REF] assume that welfare is adversely impacted by environmental damages resulting from groundwater depletion, and they empirically show there are non-negligible welfare gains from the management of two large aquifers in Southern Spain. Perreau et al. [START_REF] Perreau | Optimality versus viability in groundwater management with environmental flows[END_REF] also apply their HEM on the Western La Mancha aquifer in Spain, but they go a step further considering aquifer drainage, which contributes to environmental ecosystems. Damages thus result from a reduction in drainage and a piezometric level decline. Although those works are of high interest, they still contrast optimal management to other management regimes. We rather specifically investigate hydrological issues that arise when farmers exploit a multilayer aquifers system without considering future consequences (i.e., myopic behavior). This assumption, especially if we consider the issue of drying up natural springs, is more adapted to our empirical application to the NWSAS.

The NWSAS extends to a large desert and arid area, where surface water is nonexistent. Consequently, the population is highly dependent on this huge reserve of groundwater such that the use of this aquifer system has intensively increased over the last 40 years, mainly for irrigation purposes. Moreover, this water system is shared between Algeria, Libya, and Tunisia. Even if there is a mechanism for the joint management of the NWSAS, the three countries still withdraw water at rates that are beyond renewal rates. The intensive exploitation of this poorly managed resource leads to a series of issues, including declines in groundwater reserves, outlet depletion, and salt intrusion. Furthermore, there is no consideration for water flow toward outlets in the Algerian and Tunisian Chotts or that of usually feeding foggaras in the oases. The risk of drying up such natural springs specifically threatens population and groundwater-dependent ecosystems. 2 Therefore, it is a challenging research to examine possible trends of abstraction in the NWSAS.

The objective of this study is two-folded. First, we extend the settings of Gisser and Sanchez [START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF] to account for the hydrological and economic features of the NWSAS. To the best of our knowledge, there is no academic application to the NWSAS. 3 More precisely, we develop a stylized HEM with two interdependent superposed reservoirs and natural drainage. Therefore, we account for vertical flows between the two aquifers and horizontal flows, both depending on water stock. Water exchanges between aquifers specifically depend on the relative difference in piezometric levels, while drainage is defined as a linear function of the piezometric level by Pereau et al. [START_REF] Perreau | Optimality versus viability in groundwater management with environmental flows[END_REF]. As such, water extraction from one reservoir will impact leakage toward the other reservoir and drainage. Few studies investigate the exploitation of multiple aquifers: Zeitouni and Dinar [START_REF] Zeitouni | Mitigating negative water quality and quality externalities by joint management of adjacent aquifers[END_REF] consider water flows between reservoirs but rather concentrate on contamination of all other connected aquifers; Roumasset and Wada [START_REF] Roumasset | Ordering the extraction of renewable resources: The case of multiple aquifers[END_REF] show in a numerical simulation on the Oahu aquifer system that all resources are simultaneously used in the long run, but the "leakier" aquifer should be used first. These studies essentially concentrate on the optimal use of multiple aquifers. Contrastingly, we consider other economic hypotheses that fit with the NWSAS: both reservoirs are exploited by myopic users with a quadratic stock-dependent cost function, and the water demand is non-stationary. Farmers consequently ignore their impact on the system as a whole.

In a second step, we calibrate parameter values using data from the NWSAS over 1955-2000. We essentially proceed in two steps. We first use the available hydrologic data to calibrate the physical parameters of our dynamic model. We then estimate the economical parameters of the instantaneous water extraction functions by using extraction and piezometric data over 1955-2000. We thus end up with a calibrated dynamic system of the evolution of the piezometric levels of the two main aquifers of the NWSAS that can be used to assess the impact of an open-access regime on the evolution of the state of the NWSAS and environmental flows. We specifically investigate time trajectories of the evolution of the amount of outflows from the system, that of piezometric levels in the two different reservoirs, and the extinction of leakages. We find that the aquifer system is exploited beyond twice the total recharge in the long run. We further find that the natural outlets of the two reservoirs might be dried before 2050. We also ask the question of the sustainability of this exploitation regime by looking at the evolution of the different water withdrawals compared to the recharge.

The paper is organized as follows. In the following section, we provide a description of the NWSAS. Section 3 describes the methodology concerning the aquifer dynamics and the water extraction decisions. Data and calibration are presented in Section 4. In Section 5 we present and discuss the simulation results. Finally,Section 6 concludes the study.

Study area

The NWSAS is a large transboundary aquifer system shared by Algeria, Tunisia, and Libya. This geological basin corresponds to a multi-layer aquifer system with depths of 3000m. The aquifer formations have been grouped into two connected homogeneous and superposed reservoirs: the mostly unconfined "Complexe Terminal"(the CT), which constitutes the upper aquifer, and the confined "Continental Intercalaire" (the CI), which is the deepest and the most extended formation of the system. The CT and the CI aquifers are separated by a semi-permeable formation with low permeability, an aquitard. Fig. 1 represents the vertical structure of the NWSAS.

Recharge occurs at the outcrops of these regional-scale aquifers representing 20% of the total surface area of the CI and most of the surface area of the CT. The Atlas Mountains and the Central Sahara outcrops represent the main recharge areas for the CI and the CT, respectively. The main discharge areas of this system are represented in Fig. 2: The Algerian and Tunisian Chotts, Foggaras (traditional hydraulic systems to capture water from aquifer) mostly in Algeria, and the outlets into the Mediterranean Sea, that is, the Gulf of Gabes and the Sirte Bay. These outflows are of high importance for downstream ecosystems. Indeed, several sites in the NWSAS are recognized as Wetlands of International Importance (Ramsar Sites), for example, the Chott areas in the northern Sahara, some oases, or the Sirte Bay in Libya. More precisely, Chott El Jerid in Tunisia or Chott Merrouane in Algeria is characteristic wetlands of the northern Sahara, where surrounded oases water needs are supplied by the outflows from the CT. In Algeria, the oasis Ouled Saïd in the southern Sahara (Foggaras area, Fig. 2) is one oasis still using Foggeras to extract water from the CT, and it represents an important habitat for migratory waterbirds. While water abstraction was quite small and mostly from the CT until the early 80's, the intensive use from both reservoirs in latter years has resulted in declines in piezometric levels. 4 Extraction rates reach nearly three times the recharge rate (AbuZeid and Elrawady [START_REF] Abuzeid | North Wester Sahara Aquifer System (NWSAS) 2012 State of the Water Report, Monitoring & Evaluation for Water In north Africa (MWINA) Project[END_REF]). The total number of withdrawal points has increased from 8,800 in 2000 to 18,160 in 2008. 5 A major part of the total withdrawals is extracted for agriculture (approximately 80%).

An important socio-economic survey launched by the Sahara and Sahel Observatory (the OSS report [START_REF] Matoussi | Aspects socio-économiques de l'irrigation dans le bassin du SASS, une meilleure valorisation de l'eau pour une gestion durable du bassin[END_REF]) reveals that the NWSAS is still inefficiently exploited by farmers using private boreholes and having free access to the water network. Moreover, property rights failures are an important issue for the NWSAS exploitation. For instance, it was observed that only 80 out of 5,600 wells in Tunisia are authorized with permits.

Intensive pumping implies several consequences. The first evidence is a depletion of the groundwater reservoirs, which in turn affects pumping costs. Higher costs represent an important issue since pumping now arises in the deeper aquifer, and we observe the end of artesianism in confined sections. Declining groundwater reservoirs storage and piezometric levels also cause discharge decrease (e.g., Gonçalves et al. [START_REF] Gonçalvès | Revisiting recharge and sustainability of the North-Western Sahara aquifers[END_REF]) affecting downstream users or ecosystems,6 but also affects the vertical flows between reservoirs. As aquifer exchanges depend on the piezometric levels difference, water moves from the aquifer with the highest piezometric level to the aquifer with the lowest level. In other words, intensive use of an aquifer may drive declines in the storage of the second reservoir, which then may affect pumping decisions.

Despite the international coordination to manage the NWSAS, the growing population in such an arid area would further pressure an already stressed system. Therefore, it is of high importance to understand how short-term private decisions will affect the whole system and the potential feedback on the use of NWSAS.

Methodology

This section introduces the HEM of private groundwater extraction for the NWSAS. We first present a stylized hydrological model that summarizes the main physical features of the NWSAS. Second, we suggest an economic model underlying extraction behaviors, given the piezometric levels of both reservoirs. We finally introduce equilibrium conditions to obtain observable extraction functions that can be estimated and integrated in our HEM.

The hydrological system

Let us consider two superposed reservoirs of different extend and depths. The first layer represents the CI, which is a confined aquifer covering the largest area denoted by A I . The second layer is an unconfined reservoir of an area of A T and relatively shallow depth. This corresponds to CT. For the sake of simplicity, we normalize the bottom of the deepest reservoir to 0, and that of the upper reservoir to h. Furthermore, the top of the unconfined layer is denoted by h, which is also the land surface elevation. We denote by h i (t), i ∈ {I, T }, the piezometric level of the CI and the CT measured in meter (m) at time t. It should be noted that a piezometric level is the elevation of the free surface of the water in a well. As the CI is a confined aquifer, it is possible that h I (t) > h T (t). 7 Each reservoir is refilled by an exogenous recharge, respectively denoted by R T > 0 and R I > 0. A proportion, β i > 0 with i ∈ {I, T }, of water, naturally flows out of these two reservoirs when piezometric levels are higher than a fixed leakage threshold, respectively hI and hI for the CI and the CT. Following Darcy's Law, the leakage is proportional to the relative difference in piezometric levels and the leakage threshold if positive or zero else. The leakage at time t therefore writes:

i (h i (t)) = β i max h i (t) -hi , 0 with i ∈ {I, T } (1)
Vertical exchanges between the CI and the CT also exist. The direction of groundwater movements across the semi-permeable layer depends on the difference in the piezometric levels: water flows from the aquifer with the highest piezometric level to the lowest one. Using Darcy's law again, we even know that the amount of water moving from one aquifer to the other is a proportion, α > 0, of the difference between both piezometric levels. We can furthermore notice that this amount of water will positively affect changes in the water balance of the aquifer with the lowest piezometric level while it will adversely impact the other aquifer. Vertical flows at time t from the CT to the CI are thus given by:

e (h I (t), h T (t)) = α (h T (t) -h I (t)) (2)
Of course, the exchange from the CI to the CT is the negative of e (h I (t), h T (t)). So, if there is no pumping, a standard water balance provides the following dynamics of the two piezometric levels, h T (t) and h I (t):

s I A I ḣI (t) = R I + e (h I (t), h T (t)) -I (h I (t)) ≡ N I (h I (t), h T (t)) s T A T ḣT (t) = R T -e (h I (t), h T (t)) -T (h T (t)) ≡ N T (h I (t), h T (t)) (3)
Hydrological parameters are summarized in Table 1.

Pumping behaviors and water demand

Following Gisser and Sanchez [START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF], we assume a perfectly competitive water market where a representative water supplier sells water to any water users at p(t) the current water price. Water users do not care about the reservoir this water is extracted. Being mainly farmers, they demand a quantity of water that equates the on-farm marginal return of water to p(t). 8Symb. We even assume, like Gisser and Sanchez [START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF], that this instantaneous water demand is linear in price, but we introduce a time non-autonomous constant. This term captures the increase of the water demand across time that takes place in the NWSAS area during the last fifty years (see section 4). More specifically, our water demand is given by:

Interpretation

w D (p(t), t) = max {a(t) -bp(t), 0} (4)
where w D is the global water consumption in m 3 /s at time t, with a(t) the maximal water demand at time t and b the sensitivity of the water demand to a instantaneous price change p(t). We finally assume that after using a part γ ∈ (0, 1) of the total amount of water consumption flows back to the upper layer of the NWSAS, that is, to the CT.

The water supplier sells the total amount of water on a single market, but she nevertheless arbitrates between extractions from the CI and the CT resp. denoted by w I and w T measured in m 3 s -1 . The optimal extraction strategy follows from instantaneous profit maximization integrating cost functions specific to each aquifer. This is consistent with the study of Gisser and Sanchez [START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF] as pumping cost depends on the depth from which groundwater is extracted, and both aquifers do not share the same piezometric level. We even extend these standard pumping cost functions to consider all the other costs to supply water (e.g., conveyance costs, maintenance...). More precisely, we consider (i) pumping costs given by the energy to lift one cubic meter of water over one meter, c h -h i w i with i ∈ {I, T }, and c the energy cost, (ii) and a quadratic exploitation cost, k i 2 w 2 i . Total water production cost is therefore summarized by:

C i (h i , w i ) = c h -h i w i + k i 2 (w i ) 2 withi ∈ {I, T } (5)
As the water market is perfectly competitive, the representative water supplier is assumed to be a price-taker and chooses the amount of water pumped in the two respective reservoirs based on the following instantaneous profit maximization: max

(w I ,w T )≥0 p × i∈{I,T } w i -i∈{I,T } C i (h i , w i ) (6)
The computation of the first order conditions immediately show that extraction in each reservoir requires that the water price is equal to the marginal production cost.

p ≤ c h -h i (t) + k i w i ∀i ∈ {I, T } (= 0 when w i > 0) (7)
We get the following levels of extraction in the two respective reservoirs at time t:

w i (p, h i ) = 1 k i max p -c h -h i , 0 with i ∈ {I, T }, (8)
and the total water supply: So far, we just have characterized water demand, the respective amount of water extracted in the two reservoirs, and the total quantity of water production. All that variables will be introduced in the hydrological model.

w S (p, h I , h T ) = i∈{I,T } 1 k i max p -c h -h i , 0 (9)

Hydro-economic model of private groundwater extraction

We now integrate the hydrological model and the economic model by introducing the economic equilibrium water quantities into the water dynamics (3). This equilibrium results from an instantaneous water price adjustment, which ensures a perfect match between the demand(4) and the supply (9):

w D (p(t), t) = w S (p(t), h I , h T ) (10)
From this market clearing condition, we can derive the equilibrium price p * (t, h I , h T ) and obtain the amount of water effectively extracted in both reservoirs at time t, w * i (h I , h T , t) with i ∈ {I, T } by using Eq.( 8). However, Eq.( 9) delineates several cases depending on the two piezometric levels (h I , h T ) and the time period t at which the market equilibrium is considered. The interior solution where pumping in both reservoirs jointly occurs is specifically characterized as follows. 9 From Eq.( 8), we know that extraction in reservoir i ∈ {I, T } only occurs if the equilibrium price is larger than a lower bound p min i (h i ) = c h -h i . Hence, joint extraction requires that the equilibrium price is larger than the maximum of the two respective lower bounds, p min (h I , h T ) = max i=I,T

p min i (h i ) (11)
As water demand is decreasing in price, this situation requires that the demand at that minimal price is larger than the total supply, that is, D p min , t > S(p min , h I , h T ) > 0. In this case, the market equilibrium price is, by computation, given by:

p * (t, h I , h T ) = k I k T bk I k T +k I +k T a(t) + c( h-h I ) k I + c( h-h T ) k T (12)
if and only if

p min (h I , h T ) < p * (t, h I , h T ) (13)
Water extraction in reservoir i ∈ {I, T } at the equilibrium now becomes as follows:

w i (h I , h T , t) = 1 k i p * (t, h I , h T ) -c h -h i with i ∈ {I, T } (14)
and water consumption as the sum of the equilibrium water extraction:

w (h I , h T , t) = w I (h I , h T , t) + w T (h I , h T , t) (15)
Finally, we insert water extraction equilibria into water dynamics (see Eq.( 3)), accounting for return flows to characterize our HEM for the NWSAS fully. As we have assumed that a proportion γ ∈ (0, 1) of total consumption w flows back to the upper aquifer (the CT), the dynamics of piezometric levels of our HEM are given by:

s I A I ḣI (t) = N I (h I (t), h T (t)) -w I (h I (t), h T (t), t) s T A T ḣT (t) = N T (h I (t), h T (t)) -(1 -γ) w T (h I (t), h T (t), t) + γw I (h I (t), h T (t), t) (16)

Calibration & parameters estimation

This section calibrates the different parameters to apply the theoretical HEM to the NWSAS. We mainly use two data sources: the work of Baba Sy [START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF] on the NWSAS and the internal report of the Observatoire du Sahara et du Sahel (OSS [START_REF] Oss | Système aquifère du Sahara Septentrional, Hydrogéologie[END_REF]). We first evaluate the hydrological parameters (Table 1) and estimate economic parameters (Table 2) using observations on pumping rates over 1955-2000.

The hydrological parameters

NWSAS is one of the largest aquifer systems with a relatively small natural recharge. The CI extends over a greater area with A I = 1 million Km 2 (OSS [START_REF] Oss | Système aquifère du Sahara Septentrional, Hydrogéologie[END_REF]), whereas the CT covers an area of A T = 650, 000 Km 2 (Gonçalves et al. [START_REF] Gonçalvès | Quantifying the modern recharge of the "fossil" Sahara aquifers[END_REF]). The analysis of the groundwater budget of Baba Sy ([4]) even provides some information on the recharge. 10 We deduce a constant recharge rate of 11.4m 3 s -1 for the CI and of 18.2m 3 s -1 for the CT.

As our model is a schematic representation of the entire domain of the NWSAS, most spatial data are averaged over the entire domain. As such, we derive an average elevation of the NWSAS, h = 992m, using the topographic maps. Maps from Baba Sy [START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF] are scanned using ArcGIS software to compute a mean value of the bottom of both reservoirs and the average thickness of the aquitard in between. 11 For clarity, we choose to normalize the bottom level of the deeper aquifer (the CI) to zero, h I = 0. We then get a value of h T = 555m for the CT and the average thickness of the aquitard, e = 220m. This last observation leads to a calibration of the coefficient of vertical exchanges, α, using Eq.( 2) between the two reservoirs. Using Darcy's Law, α = KA T e 0.03m 2 s -1 where the permeability coefficient K = 10 -11 ms -1 and A T = 6.5 × 10 11 m 2 .12 . We apply the same methodology to deduce the piezometric levels in 1950 and 2000. 13 . Accounting for our normalization: these piezometric levels are given by: h I (1950) = 952.4, h T (1950) = 841.9, h I (2000) = 921.4, and, h T (2000) = 821.9 [START_REF] Esteban | The problem of sustainable groundwater management: the case of La Mancha aquifers, Spain[END_REF] The set of piezometric levels of 1950, which is often identified to the steady-state without pumping, are used as an initial condition for our piezometric dynamics ( 16), whereas the observed piezometric levels of 2000 are used as control points. The leakage coefficients, β i , and the thresholds, hi , present in the leakage equations ( 1) also are calibrated by using this reference piezometric levels of 1950 and 2000 associated to the estimates of the leakages at the same years provided by Baba Sy's water balances [START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF] of 1950 and 2000. 14 They are given by: I (1950) = 8m 3 s -1 , T (1950) = 13.7m 3 s -1 and I (2000) = 5.4m 3 s -1 , T (2000) = 4.5m 3 s -1 . By linear interpolation, we get hI , hT (857.02m, 812.12m) and (β I , β T ) (0.084m 2 s -1 , 0.46m 2 s -1 ).

Finally, we estimate the two storativity coefficients, s I and s T , using the above estimated parameters on leakage, leakance coefficient, data on the piezometric maps for 1950 and 2000, and on water pumping available in the internal report of the OSS ( [START_REF] Oss | Système aquifère du Sahara Septentrional, Hydrogéologie[END_REF]) 15 and reproduced in Table 3. These data make it possible to reconstitute the temporal evolution of the piezometry between 1950 and 2000 by using successive water balances and adjusting the storativity coefficients by minimizing the square of the difference, in 2000, between observed and simu-lated piezometric levels. Using Excel's solver, we get s I = 0.00051 and s T = 0.00073 as well as the time series of the piezometric levels given in Table 3.

Years w I in m 3 s -1 w T in m 3 s -1 w in m 2 s -1 h I in m h I in

The economic parameters

We now turn to the calibration of the parameters of water demand (4) and exploitation costs [START_REF] Bekchanov | Systematic review of Water-Economy modeling applications[END_REF]. First, following Gonçalves et al. [START_REF] Gonçalvès | Quantifying the modern recharge of the "fossil" Sahara aquifers[END_REF], we assume that 15% of the pumping returns to the aquifer system, that is, γ = 0.15, and we make the assumption that almost all return flow reaches the CT. We secondly compute the cost of pumping water, c, using the formula to assess the energy in Joule (J) requirement to lift one cubic meter of water over one meter: ρg e , with ρ = 10 3 kgm -3 the volumetric mass of one m 3 of water, g 9.81ms -2 the gravity acceleration and e the pump average efficiency fixed at 75%. If we now consider that a liter of diesel oil represents an energy value equal to 1 3.787.10 7 J and that the average price of diesel oil is p 0.23e per liter, 16 We can easily assess the pumping cost c 7.944.10 -5 e(m m 3 ) -1 . It also remains to specify the non-stationary part of the water demand, a(t). This allows capturing the increase of the water demand observed in the NWSAS area during 1950-2000. The total groundwater consumption over this period (see column 3 in Table 3 or Fig. 3) largely suggests the use of a logistic function given by a(t) = a 0 + a 1 1+a 2 e -rt . Consequently, four new parameters (a 0 , a 1 , a 2 , r) need to be estimated in addition to the other parameters (b, k I , k T ). All these parameters are needed to characterize the pumping levels, w i (h I , h T , t) i ∈ {I, T ) (see Eq.( 14)), which also depend on the piezometric levels of the CI and the CT reservoirs. We will use the set of data in Table 3 to estimate these parameters by solving a nonlinear multi-curve fitting problem based on least-square minimization. As we only have a few data and we know that this estimation method is sensitive to initial conditions, we run a set of 500 estimates using Matlab with randomly chosen initial conditions, and we select the set of parameter values that provides the best curve fit outlined in Fig. 4:

Table 5 summaries the estimation of the different economic parameters.

Results

The calibration of the parameters is now used to simulate the evolution of the NWSAS after the period 1950-2000. We proceed in two steps. We first discuss changes in the piezometric levels of the CI and the CT with the subsequent effects on the natural discharge of both aquifers and vertical fluxes. We then address the question of the sustainable exploitation of this aquifer system.

On the evolution of the hydrological system

We first characterize the time trajectories of the piezometric levels of both reservoirs, the CI and the CT. Plugging parameter values of Tables 4-5 in the dynamic system [START_REF] Duarte | Waer use in Spanish economy: An input output approach[END_REF], we get the following differential system: 10

The first linear part of this equation mainly reflects the effect of the recharge, the leakance, and the pumping activity (net of the return flow) on the piezometric levels. The nonautonomous term captures the exogenous evolution of the water demand during the period, while the last non-linear term reflects the potential regime shifts induced by the interruption of the two aquifers' natural leakage. By setting the initial conditions to the almost natural piezometric levels observed in 1950 (see conditions ( 17)), we find a numerical solution to this system using Matlab, which are depicted in Fig. 5. We easily observe a significant decline in the two piezometric levels after the mid-1970s. This is consistent with the intensive use of the NWSAS, which started in the same period (see Fig. 4). Such an increase in the exploitation of the system has caused a significant drop in the piezometric levels over the observation period (1950-2000: The CI lost approximately 30 meters and the CT 20 meters. 1.8526 -0.13075 -0.13075 2.9090

h I (t) h T (t) (19)
Such exploitation will lead to a drawdown of approximately 150 meters and 100 meters of the piezometric levels of resp., the CI, and the CT, over 1950-2100. Moreover, no steady state appears to be reachable. Furthermore, we even observe in Fig. 5 that the piezometric levels fall below the CI and CT leakage thresholds. Our model, based on average piezometric levels, suggests that the natural discharges of the CT will dry up after 2010, a situation we locally observed with depleted wells. The discharge of the CI should also be completely dried up in 2045, as we can notice in Fig. 6. As soon as the CT reservoir became intensively used, natural leakage sharply declined after a period of the increase until the mid-1970s. This thus suggests that intensive pumping has caused this issue. This is less obvious for the CI. The leakage has diminished since 1950, while this reservoir has been exploited later than the CT. This primary results from a natural process. Indeed, we easily observe that the piezometric level of the CI is always higher than that of the CT. As the vertical fluxes are driven by the differential in the piezometric levels, the CI always feeds the CT. As such, the piezometric level of the CI declines and the natural discharge. As pumping rates have started increasing in the early 1980s, we nevertheless observe a higher decline in the natural drainage of the CI. Even if the decline is slower, natural discharges from the CI should also be dried up. Thus, we can assert that all the remaining natural outlets of the NWSAS are now at risk, with potentially irreversible damage for the downstream users, especially the Saharan ecosystems that benefit from this water.

Declines in piezometric levels finally alter vertical fluxes. Fig. 7 shows that leakance from the CI toward the CT has always decreased. We can distinguish three phases. Over 1950-1990, we observe a steady decline since the differential in piezometric levels decrease. This is not so surprising because the piezometric level of the CI declines faster than the level of the CT (3). Indeed, a higher level of water flows out of the reservoir [START_REF] Berrittela | The economic impact of restricted water supply: A computable general equilibrium analysis[END_REF]. The second phase is characterized by a sharp drop in leakance. This more or less corresponds to the period of the intensive use of the CI and should finish when there is no natural drainage anymore. During this period, the piezometric level of the CI should decline more rapidly than that of the CT, leading to decrease in the differential of both levels, and consequently the level of vertical exchanges of water. The third and last phase is characterized by a stable level of leakance. This thus means that both piezometric levels are decreasing at the same rate.

This last comment calls for greater attention. Indeed, the time evolution of the two piezometric levels results from different levels of inflows and outflows. Typically, the less naturally replenished reservior (the CI) feeds the reservoir that receives the higher level of natural recharge (the CT). Moreover, the CT reservoir receives part of irrigation water, that is, another amount of water extracted from the CI. A perusal analysis of the water balances should thus provide interesting information on the degree of over-exploitation of each reservoir.

Evolution of the water balances

We now turn to the analysis of the water balance of the NWSAS and that of each reservoir. This allows assessing whether an aquifer is over-exploited by contrasting aquifer inflows with outflows. At the global scale, we will thus compute all water inflows, that is, the sum of pumping net of return flows and natural discharges, to compare with the sum of the two natural recharges. Fig. 8 shows that the exploitation of the NWSAS is unsustainable. Total pumping, which was in the mid-1070s far below the total recharge of the NWSAS, increased in the late 1990s to twice the recharge so that the total water outflows peak at 2.5 times a recharge. Today, total pumping appears to reach a steady-state at 62m 3 s -1 and total water outflows even decrease due to the gradual disappearance of natural discharges. However, this potential stabilization of water withdrawals is still largely unsustainable as the total recharge is 29.6m 3 s -1 .

However, the global analysis hides the heterogeneous situation of both reservoirs. Indeed, water imbalances are larger in the CI than in the CT. This mainly comes from leakance and return flows. Fig. 9 shows that an amount equivalent of three times the recharge runs out the reservoir at the end of the nineties. We easily observe that pumping has drastically increased from 5m 3 s -1 in the mid-1970s to an, a priori, stable level at 28m 3 s -1 , nowadays. We even know that a minimum amount of 2m 3 s -1 flows to the CT reservoir, and the natural discharge only stops in 2045. Given that the natural recharge is approximately 11.8m 3 s -1 , we easily understand that the stock of the CI is tapped. More precisely, the peak water outflows arose approximately 2000 and were well above the triple of the recharge. Even when all springs are dried up by about 2050, the volume of outflows (pumping and leakance) will continue to be more than two and half times the recharge. The situation of the CT is less obvious because the reservoir is refilled by return flows and drainage from the CI. Water imbalances are slightly acute. Fig. 10 summarizes all the outflows in the positive quadrant and all the inflows in the negative quadrant (except the natural recharge). We can first observe that the pumping level was multiplied by a factor of 4.5 from the mid-1970s to the present. We also know that there is no natural drainage anymore. This improves the balance at, of course, a certain environmental cost. Meanwhile, inflows always come from the CI because of the positive differential in the piezometric levels.

As the CT is the upper reservoir, the majority of return flows reach the reservoir. In addition to the natural recharge of 18.2m 3 s -1 , the CT will collect in the next 40 years a water inflow, on average, of 13.5m 3 s -1 . This reduces the pressure on water stock. Nevertheless, total outflows remain larger than twice the recharge.

Concluding remarks

This study contributes to the literature on applications of HEMs. We studied the state of groundwater resources in the NWSAS, a multilayer aquifer system, which has been intensively exploited for several years. We developed a theoretical model to represent specific features of that groundwater basin and provided a numeric simulation based on a calibration of the different parameters. To construct this model, we borrowed several elements of the hydrologic analysis of the NWSAS (e.g., Baba Sy [START_REF] Sy | Recharge et paléorecharge du Système Aquifère du Sahara Septentrional[END_REF], Gonçalvès et al [START_REF] Gonçalvès | Quantifying the modern recharge of the "fossil" Sahara aquifers[END_REF], [START_REF] Gonçalvès | Revisiting recharge and sustainability of the North-Western Sahara aquifers[END_REF]. Based on these results, we first specified water dynamics for each reservoir, but we accounted for the existence of vertical water flows between layers. In a second step, inspired by the economic literature that tests the robustness of the Gisser-Sanchez effect (e.g., Brill and Burness [START_REF] Brill | Planning versus Competitive Rates of Groundwater Pumping[END_REF], Tomini [START_REF] Tomini | Is the Gisser and Sánchez model too simple to discuss the economic relevance of groundwater management?[END_REF]), we adopted a growing demand to capture the development of agricultural activity in the region and population growth. Moreover, we introduced a more general cost function to account for pumping costs and other costs associated with water delivery. In a numerical simulation of this integrated model, we highlighted sharp declines in piezometric levels, as well as natural outflows. We especially showed that all outflows will be dried up in a short horizon. We also showed that the actual groundwater exploitation regime is far to be sustainable in the long run.

Even if we do not characterize the optimal use of such a system, we observe that intensive exploitation causes several environmental issues, especially a change in the resource regime since leakages disappear in short term as piezometric levels drop. This thus causes stress on groundwater-dependent ecosystems, local populations, and agriculture. Such observations associated with an unsustainable exploitation regime call for public policies to adequately use water from the ground and protect natural outflows.

Nevertheless, we did characterize the optimal policy to measure benefits from intervention the impacts induced by other property rights regimes. Therefore, it would be interesting to extend this analysis to design socially-optimal exploitation of the NWSAS. This will allow us to examine whether natural outflows will be dried up, especially if the social planner considers damages on ecosystems. As the resource is shared in three countries, it would also be interesting to consider strategic externality and wonder how to design cooperation to help the ongoing cooperative process. Finally, we consider a schematic representation of the entire domain. We could also consider spatial heterogeneity as some areas have only access to one reservoir.

Figure 1 :

 1 Figure 1: Schematic cross-section of the NWSAS

Figure 2 :

 2 Figure 2: Extension of the CT and CI aquifers with major outflow areas

Figure 3 :

 3 Figure 3: Total Groundwater consumption over 1950-2000

Figure 4 :

 4 Figure 4: Observed and estimated water extraction paths

Figure 5 :

 5 Figure 5: Evolution of the piezometric levels in the NWSAS

Figure 6 :

 6 Figure 6: Evolution of the leakages

Figure 7 :

 7 Figure 7: Evolution of the leakance CI->CT

Figure 8 :

 8 Figure 8: Global water balance of the NWSAS

Figure 9 :

 9 Figure 9: Water balance of the CI

Figure 10 :

 10 Figure 10: Water balance of the CT

Table 1 :

 1 Hydrological Parameters

	Units Symb.	Interpretation	Units

Table 2

 2 summarizes all economic variables and parameters.

		Economic Variables			Economic Parameters
	Symb.	Interpretation	Units	Symb.	Interpretation	Units
	w(t)	Water consumption	m 3 s -1	a(t)	Maximal demand at t	m 3 s -1
	w I (t) Extraction from the CI	m 3 s -1	b	Slope of the demand	(m 3 ) 2 e -1 s -1
	w T (t) Extraction from the CT	m 3 s -1	c	Coef. of the lift cost	e(m m 3 ) -1
	p(t)	Water price at time t	e(m 3 ) -1	k I	Coef. of the quad. cost (CI)	e(m 3 ) -2
				k T	Coef. of the quad. cost (CT)	e(m 3 ) -2
				γ	Rate of return flow	-

Table 2 :

 2 The economic variables and parameters

Table 3 :

 3 Pumping rates, water consumption and piezometric levels over

	m

Table 4

 4 summarizes all the parameter values.

	Symb.	Description	Values for the CI Values for the CT
	A i h	Area Surface elevation	10 12 m 2 992 m	6.5 × 10 11 m 2 992 m
	h	Bottom of the reservoir	0 m	555 m
	s i	Storativity	0.00051	0.00073
	R i	Recharge	11.4 m 3 s -1	18.2 m 3 s -1
	hi	Leakage threshold	857.02 m	812.12 m
	β i	Leakage coef.	0.084 m 2 s -1	0.46 m 2 s -1
	α	Leakance coef.	0.033 m 2 s -1

Table 4 :

 4 Hydrological parameter values

Table 5 :

 5 Estimates of Economic parameters

	Symb.	Values	Symb.	Description	Values
	Drift of the demand	b	Slope of the demand	51.721(m 3 ) 2 e -1 s -1
	a 0	116.445 m 3 s -1	c	Coef. of the lift cost	7.944.10 -5 em -1 m -3
	a 1	557.19 m 3 s -1	k I	Coef. of the quad. cost (CI)	0.4108 e(m 3 ) -2
	a 2	1176.8	k T	Coef. of the quad. cost (CT)	0.2546 e(m 3 ) -2
	r	0.241 t -1	γ	Rate of return flow	0.15

Globally speaking, optimization-based models and simulation-based models are often used to analyze water issues at a regional or hydrological basin scale, whereas economy-wide models are used to examine the impacts of water use on a larger scope of the economy (Berrittella et al.[START_REF] Berrittela | The economic impact of restricted water supply: A computable general equilibrium analysis[END_REF], Dinar[START_REF] Dinar | Economy-wide implications of direct and indirect policy interventions in the water sector: Lessons from recent work and future research needs[END_REF], or Duarte and Cholitz[START_REF] Duarte | Waer use in Spanish economy: An input output approach[END_REF], or Llop[START_REF] Llop | Economic impact of alternative water policy scenarios in the Spanish production system: An input-output analysis[END_REF]).

There are several different groundwater-dependent ecosystems such as oases, sabkhas or Chotts. We also find Ramsar sites in the NWSAS.

Several studies applied to the NWSAS are developed by the Sahara and Sahal Observatory in the context of two international projects:the SASS project (2010-2015) (see http://www.oss-online.org/en/north-western-sahara-aquifer-system-sass) and the NEXUS projet 2016-2019 (see http://www.oss-online.org/en/acting-cooperation-water-sector-mediterranean-nexus).

Total extractions increased from about 1 km3 in 1980 to 2.5 km3 in 2000 (Baba Sy[START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF], OSS[START_REF] Oss | Système aquifère du Sahara Septentrional, Hydrogéologie[END_REF])

We exactly know that the repartition was in 2000 as follows: 3,500 points in the CI, and 5,300 in the CT.

Many outflows are under stress, including foggaras in Algeria, the outlet discharges in the Gulf of Gabes in Tunisia, or the Chotts in Algeria and Tunisia.

This corresponds exactly to the situation observed in the NWSAS. It is even possible that h I (t) is higher than the ground level (i.e., an artesian behavior), but this last case never occurs under our calibration using bassin-scale average piezometric levels and homogeneous reservoir description.

The reader may wonder whether the existence of a water market is relevant to our case study in which farmers often own their wells. However, as long as this market is competitive, myopic water allocation strategies can be separated from myopic water extraction decisions. This observation keeps our model close to the one of Gisser and Sanchez[START_REF] Gisser | Competition versus Optimal Control in Groundwater Pumping[END_REF] and allows comparisons.

Other cases may exist where pumping occurs in a single reservoir. However, the joint use of both reservoirs fit with our application to the NWSAS. Moreover, the condition (13) is satisfied for all t in the simulation. We therefore choose to omit these additional cases from our study.

For details, see Table10-5 p.157 in Baba Sy[START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF].

More precisely, we use the maps of top of the CI and the CT (resp. Fig.5-21 p.77 and Fig.5-24 p.78) and maps of the thickness of the CI and the CT (resp. Fig.3-7 p.49 and Fig.3-9 p.51) of Baba Sy [4].

Baba-Sy[START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF] assesses that the permeability coefficient varies between 10 -9 to 10 -12 ms -1 , we arbitrarily choose an intermediate value.

We specifically use maps from Baba Sy (Maps 8-1 p.118 and 8-2 p.119) to compute the piezometric levels in 1950 and maps 12-13 p.117-118 in OSS for 2000([37]).

For details see Table10-5 p.157 in Baba Sy[START_REF] Baba-Sy | Holocen and present recharge of the Saharan aquifers: A numerical modelling study[END_REF].

 (p.103-104) and Tables

16-18 (p.106-108), respectively, provide information on extraction rates in the CT and the CI at country levels. We thus aggregate data to get total pumping rates over the period1955-2000.

This diesel oil price corresponds to an average price in Algeria, Tunisia, and Libyan over 1991-2016, according to the world bank data. This choice seems rather arbitrary, but the pumping cost coef. is so small that changes in the diesel oil price do not change the results.

Hubert Stahn and Agnes Tomini acknowledge the ANR financial support (Grant ANR-17-EURE-0020), the Excellence Initiative of Aix-Marseille University -"A*Midex" and the ANR project GREEN-Econ (Grant ANR-16-CE03-0005).