# AQUACOAST: A Simulation Tool to Explore Coastal Groundwater and Irrigation Farming Interactions

Authors: Jaime Martínez-Valderrama; Javier Ibáñez; Francisco J. Alcalá
Journal: Scientific Programming
 (2020)

---

## Abstract

In the framework of coastal groundwater-dependent irrigation agriculture, modelling becomes indispensable to know how this renewable resource responds to complex (usually not conceptualized nor monitored) biophysical, social, and economic interactions. Friendly user interfaces are essential to involve nonmodeling experts in exploiting and improving models. Decision support systems (DSS) are software systems that integrate models, databases, or other decision aids and package them in a way that decision makers can use. This paper addresses these two issues: firstly with the implementation of a System Dynamics (SD) model in Vensim software that considers the integration of hydrological, agronomic, and economic drivers and secondly with the design of a Venapp, push-button interfaces that allow users access to a Vensim model without going through the Vensim modelling environment. The prototype designed, the AQUACOAST tool, gives an idea of the possibilities of this type of models to identify and analyze the impact of apparently unrelated factors such as the prices of cultivated products, subsidies or exploitation costs on the advance of saltwater intrusion, and the great threat to coastal groundwater-dependent irrigation agriculture systems.


---

## Body

## 1. Introduction

A large proportion of the world’s population depends on groundwater resources [1]. Groundwater accounts for as much as 33% of the global freshwater demand and is increasing [2]. This resource is especially crucial in regions where the surface water network does not cover a fraction of the freshwater demand at least, such as in many coastal drylands [3]. In these areas, aquifers play a critical role in sustaining the economy and the environment [4], but water scarcity has led to high groundwater abstraction rates to supply the increasing urban, tourism, industrial, and agriculture demands, thus adding stress to groundwater bodies [5, 6]. The cases with alarming signs of groundwater quantity depletion and quality degradation do not stop increasing [7].Groundwater degradation in coastal drylands is often the result of complex interactions between global driving forces and human-induced actions [8, 9]. Global climate scenarios forecast high evaporation rates and lower high-salinity aquifer recharge rates [10]. Some human-induced actions may contribute to this initial “climatic” groundwater salinity, such as pumping inducing saltwater intrusion and mobilisation of brines [4, 11], high-efficiency irrigation systems producing high-salinity, low irrigation return rates [12], and misuse of fertilizers and agrochemicals, among others. For instance, the resulting groundwater salinity in many southern Europe coastal drylands exceed the standards of quality that the European Drinking Water Directive [13] recommends for human consumption and health, crop production, and survival of the ecosystems biodiversity [12].In this context, agriculture has invariably had the highest freshwater demand [14, 15], contributing at growing of local societies, but accounting for 70% of withdrawals [16, 17]. The exposition of coastal aquifers to saltwater intrusion (SWI) after disruption of natural conditions [18, 19] had reduced crop surface and therefore profit in many regions, with dramatic consequences for food production and social inequality. Users often forget that seawater is another water balance component. An unplanned exploitation is the cause behind the alteration of the fragile freshwater-saltwater equilibrium with typical consequences such as freshwater-saltwater mixing zone displacement inland and reduction of crop production. These anticipated negative indicators may end with abandonment of land devoted to irrigated crops.Understanding coastal groundwater and irrigation farming interactions is a challenge because the equilibrium between sustainable pumping and agricultural profit is subjected to a fragile threshold determined by the crop-limiting salinity. This paper aims to contribute to this topic by means of a user-friendly programming tool that combines hydrological, agronomic, and economic components of the system [20, 21]. This multidisciplinary tool must deal with the complex network of variables that conform feedbacks loops subjected to delays and nonlinear relationships. For this, SD modelling is a very suitable approach. When the scope of the study goes beyond the academic field and aims to involve stakeholders and land-use planners, it is essential to link these models to DSS. These are software systems that integrate models, databases, or other decision aids and package them in a way that decision makers can use [22].The perception of the complex interaction between the sociocultural, economic, and biophysical components of the systems is the cause of the growing interest in environmental DSSs [22–24]. This paper presents a Vensim application (Venapp) for the AQUACOAST simulation model, which aims to evaluate the sustainability of coastal groundwater-dependent agriculture. As described, the great threat to coastal aquifers is saltwater intrusion triggered by excessive freshwater pumping. Analyzing under what scenarios such a disaster may occur is the purpose of this simulation tool.


## 2. Materials and Methods

### 2.1. System Dynamics Modelling

SD modelling [25] is an ideal framework to integrate variables coming from different disciplines because it is based on a holistic view of the system under study and considers the model structure as a whole. SD modelling embodies system thinking through the implementation of differential equation systems. It is grounded in the theory of nonlinear dynamics and feedback control developed in mathematics, physics, and engineering.A SD model is a stock-and-flow structure of first-order differential equations, which is commonly deterministic and used to generate the time trajectories of the model variables under different simulation scenarios. The model structure is made up of causal feedback loops, which include nonlinear relationships and delays, and constitute a holistic and easily overlooked because of its behavior [26].Strong points of SD may be summarized as follows [27]: (i) SD models are useful learning tools that help improve system understanding and foster system thinking skills and knowledge integration for modelers and end users. This is true even for their first stage of development, i.e., as conceptual models. (ii) Advances in the development of high-level dynamic modelling software platforms have made computational SD modelling more accessible. (iii) SD literature has made rich contributions to approaches informing on the modelling process, including data collection methods, knowledge elicitation/mapping techniques, and policy analysis.The integrated models that allow the implementation of SD are essential to address environmental problems [28–30]. The degradation of the environment is forged by the imbalance created between the natural rate of regeneration of resources and their rate of exploitation, which is due to socioeconomic factors. In order to understand how economic systems are decoupled from natural ones and to seek sustainable solutions, it is essential to analyze the dynamic interactions between economic and natural biophysical systems [21, 31].Advances in the development of high level dynamic modelling software platforms, such as iThink (isee systems,http://www.iseesystems.com), Vensim (Ventana Systems, http://www.vensim.com), and Powersim Studio (Powersim Software AS, http://www.powersim.com), have rendered computer SD models very accessible (even with minimal technical background). These applications are often designed as communication layers: user interface, stock-flow, mathematical equations, and programming code [27].The important feature of these modelling tools is that they allow the user to write a model, even a complex one, without writing a single equation. The advantage is that the model diagram is focused on the relationships among the different variables and processes in the system and can be understood by a person with a good background in the problem domain of the model, but with little knowledge of mathematics. Models are designed “visually” by connecting blocks and aggregating them at different levels of resolution [22] what makes these tools very useful in the model prototyping phase and for educational purposes.From all available options, we chose Vensim [32]. Although the software lacks certain possibilities, for example, does not support spatial treatment and its graphical interface is improvable [33], it is very suitable for the purpose we aim, to understand the typical functioning complexity of a coastal groundwater-dependent agricultural system subjected to different water-competing sectors and land-use policies. Complexity departs from the lack of data in many groundwater systems [34], thus complicating numerical modelling. There may be unknown essential information, such as the actual groundwater stored in an aquifer. In addition, cause and effect are separated in time and space which adds difficulty in conceptualizing basic hydrological processes. For instance, aquifer recharge may take place tens to hundreds of kilometers from the discharge area and may occur months to years before this water is pumped [4, 8].


### 2.2. Vensim Applications: Venapps

Interaction with simulation models has improved considerably in recent years. However, it can still be difficult for a nonexpert user to deal directly with programs such as Vensim. Friendly user interfaces are essential to involve nonmodeling experts in exploiting and improving models. The VensimDSS® software [35] includes the possibility of developing simple Venapps to visualize the model structure, examine its causality, simulate the model behavior, and visualize its results [36–38].Venapps are simplified, push-button interfaces that allow users access to a Vensim model without going through the Vensim modelling environment. A Venapp uses a model and a set of rules for interacting with the model to give users simplified access to that model. To the user, the Venapp appears as a series of buttons, menus, or a sequence of screens allowing him or her to use and analyze the model in straightforward and meaningful way.The goal of these interactive simulators is to provide the user with the following tasks: examine the model structure, establish simulation scenarios, simulate the model, and analyze the impact of decisions on the behavior of the system under study. The VensimDSS environment includes an extension called Vensim application (VenappTM) which allows the creation of graphical user interfaces in the file format “vpa” (Vensim packaged application). Both model and Venapp can be run with the software “Vensim Model Reader” which can be downloaded from the Vensim® homepage for free (http://www.vensim.com/freedownload.html).


## 2.1. System Dynamics Modelling

SD modelling [25] is an ideal framework to integrate variables coming from different disciplines because it is based on a holistic view of the system under study and considers the model structure as a whole. SD modelling embodies system thinking through the implementation of differential equation systems. It is grounded in the theory of nonlinear dynamics and feedback control developed in mathematics, physics, and engineering.A SD model is a stock-and-flow structure of first-order differential equations, which is commonly deterministic and used to generate the time trajectories of the model variables under different simulation scenarios. The model structure is made up of causal feedback loops, which include nonlinear relationships and delays, and constitute a holistic and easily overlooked because of its behavior [26].Strong points of SD may be summarized as follows [27]: (i) SD models are useful learning tools that help improve system understanding and foster system thinking skills and knowledge integration for modelers and end users. This is true even for their first stage of development, i.e., as conceptual models. (ii) Advances in the development of high-level dynamic modelling software platforms have made computational SD modelling more accessible. (iii) SD literature has made rich contributions to approaches informing on the modelling process, including data collection methods, knowledge elicitation/mapping techniques, and policy analysis.The integrated models that allow the implementation of SD are essential to address environmental problems [28–30]. The degradation of the environment is forged by the imbalance created between the natural rate of regeneration of resources and their rate of exploitation, which is due to socioeconomic factors. In order to understand how economic systems are decoupled from natural ones and to seek sustainable solutions, it is essential to analyze the dynamic interactions between economic and natural biophysical systems [21, 31].Advances in the development of high level dynamic modelling software platforms, such as iThink (isee systems,http://www.iseesystems.com), Vensim (Ventana Systems, http://www.vensim.com), and Powersim Studio (Powersim Software AS, http://www.powersim.com), have rendered computer SD models very accessible (even with minimal technical background). These applications are often designed as communication layers: user interface, stock-flow, mathematical equations, and programming code [27].The important feature of these modelling tools is that they allow the user to write a model, even a complex one, without writing a single equation. The advantage is that the model diagram is focused on the relationships among the different variables and processes in the system and can be understood by a person with a good background in the problem domain of the model, but with little knowledge of mathematics. Models are designed “visually” by connecting blocks and aggregating them at different levels of resolution [22] what makes these tools very useful in the model prototyping phase and for educational purposes.From all available options, we chose Vensim [32]. Although the software lacks certain possibilities, for example, does not support spatial treatment and its graphical interface is improvable [33], it is very suitable for the purpose we aim, to understand the typical functioning complexity of a coastal groundwater-dependent agricultural system subjected to different water-competing sectors and land-use policies. Complexity departs from the lack of data in many groundwater systems [34], thus complicating numerical modelling. There may be unknown essential information, such as the actual groundwater stored in an aquifer. In addition, cause and effect are separated in time and space which adds difficulty in conceptualizing basic hydrological processes. For instance, aquifer recharge may take place tens to hundreds of kilometers from the discharge area and may occur months to years before this water is pumped [4, 8].


## 2.2. Vensim Applications: Venapps

Interaction with simulation models has improved considerably in recent years. However, it can still be difficult for a nonexpert user to deal directly with programs such as Vensim. Friendly user interfaces are essential to involve nonmodeling experts in exploiting and improving models. The VensimDSS® software [35] includes the possibility of developing simple Venapps to visualize the model structure, examine its causality, simulate the model behavior, and visualize its results [36–38].Venapps are simplified, push-button interfaces that allow users access to a Vensim model without going through the Vensim modelling environment. A Venapp uses a model and a set of rules for interacting with the model to give users simplified access to that model. To the user, the Venapp appears as a series of buttons, menus, or a sequence of screens allowing him or her to use and analyze the model in straightforward and meaningful way.The goal of these interactive simulators is to provide the user with the following tasks: examine the model structure, establish simulation scenarios, simulate the model, and analyze the impact of decisions on the behavior of the system under study. The VensimDSS environment includes an extension called Vensim application (VenappTM) which allows the creation of graphical user interfaces in the file format “vpa” (Vensim packaged application). Both model and Venapp can be run with the software “Vensim Model Reader” which can be downloaded from the Vensim® homepage for free (http://www.vensim.com/freedownload.html).


## 3. AQUACOAST Model

### 3.1. Overall Description

The AQUACOAST model implements the two-dimensional saltwater-freshwater interface (SFI) analytical solution proposed in [39] for an idealized unconfined (i.e., its upper boundary is the fluctuating water table) coastal aquifer, which is horizontal, enough thick, and overlying an impervious bedrock, isotropic, and hydraulically homogeneous with steady inflow (aquifer recharge) and outflow (aquifer discharge) rates (Figure 1). The variables and parameters of the AQUACOAST model concerning irrigation practices, groundwater hydrology, and groundwater allotment for agriculture have been added to the schematization of Verruijt [39]. In the steady schema of Verruijt [39], the SFI position is reached automatically depending on the fresh groundwater flow. AQUACOAST modifies this hypothesis by introducing delays in order to reach the equilibrium described by Verruijt, which magnitude differs, which differ for the intrusion phase (minor delay) and the freshening phase (longer delay) [40]. Groundwater allotment for agricultural use is the main flow component. The model focuses on determining critical distances of agronomic interest from the shoreline, such as the current and equilibrium freshwater and salinization boundaries or the irrigation farming abandonment boundary (Figure 1).Figure 1
The AQUACOAST model added to Verruijt’s schematization for an idealized unconfined coastal aquifer.hF is the height of the free groundwater surface (FGS: water table) above sea level (it is not a model parameter), hS is the depth of the saltwater-freshwater interface (SFI) below sea level (it is not a model parameter), hght is the average height of land surface above sea level, vtmz is the vertical thickness of the mixing zone, dpwl is the depth of pumping wells below sea level, CDSB is the current distance of the SB from the shoreline, CDAB is the current distance of the AB from the shoreline, CDFB is the current distance of the FB from the shoreline, EDSB is the equilibrium distance of the SB from the shoreline, EDFB is the equilibrium distance of the FB from the shoreline, lgth is the coastal aquifer length, slfw is the freshwater salinity, S is the mixing water salinity along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB), slmz is the actual average salinity across a line segment connecting FB and AB (CDAB ≤ x ≤ CDFB), slcl is the crop-limiting salinity, and slsw is the saltwater salinity.AQUACOAST is a lumped spatial model [27], and thus, its variables represent either totals or averages over the entire area modelled. For the sake of illustration, the simulation period is 50 years, an appropriate time horizon when studying environmental degradation that can only be detected when observed the evolution of “slow variables” [41]. Variables and parameters use the SI units, capital letters being variables, whereas small letters denote parameters. Tables 1 and 2 include the notation, units, and adopted values for variables and parameters, respectively, after data from [4, 8, 20, 21, 40, 42, 43] and references therein. The AQUACOAST model uses 26 straightforward formulations divided into seven sections as follows.Table 1
Notation for variables used.
VariableEquationDefinitionUnits1ARIR15, 26Current irrigated aream−2·day−1ARPI11, 13Potential irrigated aream3·ha−1·year−1ARTI13, 15Target for the irrigated area€·year−1CDAB8, 11, 14, 19Current distance of the AB from the shorelinem3·ha−1·year−1CDFB5, 7, 8, 19Current distance of the FB from the shoreline€·ha−1·year−1CDSB3, 6, 7, 8, 14Current distance of the SB from the shoreline€·m−3CHFW10, 20Current average height of the FGS above sea levelmCSTW20, 21, 24Water cost per cubic metermEDFB4, 5Equilibrium distance of the FB from the shorelinem3·ha−1·year−1EDSB2, 3, 6, 14Equilibrium distance of the SB from the shorelinem3·ha−1·year−1EHFW9, 10Equilibrium average height of the FGS above sea levelkmEXPF13, 22Expected profit per irrigated farmkmEXWX23, 25Expected water requirement for profit maximizationkmIRDS17, 21, 23, 26Irrigation dose from groundwaterkmPFTH21, 22Profit per hectarekmPPIR16, 17, 23Contribution of direct precipitation to irrigationyear−1PUMP1, 26Groundwater pumpingdmnlRAIA14, 15Rate of abandonment of irrigated areahaS(x)7, 18Salinity of mixing water along thex-axis distance from the shoreline, CDSB ≤ x ≤ CDFBkg·m−3SLFT17, 18, 19, 24Salinity factorhaTDFW1, 2, 4, 9Net aquifer dischargehaTMBE3, 5, 6, 14Time for the SB and FB to achieve equilibriumyearWRXP24, 25Water requirement for profit maximizationmm·year−1YLDH17, 21Crop yield per surface unitkg·ha−1·year−11 – dmnl = dimensionless; ha = hectare; € = EU Euros.Table 2
Notation for parameters used.
ParameterEquationDefinitionUnits1Value2aaet1Annual actual evapotranspirationmm·year−1200appt1, 16Annual precipitationmm·year−1300arfr15, 22Area of one farmha2.84cdap7, 8, 19Calibration parameter to represent nonlinear variations in salinitydmnl3cfpm2, 4, 9Hydraulic conductivitym·day−110cfrf1Return flow coefficientdmnl0.2cfro1, 16Runoff coefficientdmnl0.2cwot20Costs of water other than energy€·m−31.4dpwl2, 4Depth of pumping wells below sea levelm100dtcp2, 4, 9Density contrast parameterdmnl0.025dtfw2, 4, 9Freshwater densitykg·m−31.000dtsw2, 4, 9Saltwater densitykg·m−31.025efir21, 24, 26Irrigation system efficiencydmnl0.8egcm20Energy required to pump one cubic meter one meterkwh·m−40.017expf22Initial expected profit per irrigated farm€·year−123,000hght20Average height of land surface above sea levelm100idqt23Irrigation quotam3·ha−1·year−120,000lgth1, 2, 4, 9, 11, 14, 19, 26Coastal aquifer lengthkm10lttf1Lateral transference from other aquifershm3·year−17.9oaag16If 1, then open-air agriculture else under plastic agriculturedmnl0ocfm13Average opportunity cost for a farmer€·yr−123,000ocfv13Coefficient of variation of the opportunity cost for a farmerdmnl0.29octh21Other costs per hectare€·ha−1·year−140,053prcp21, 24Crop price€·kg−10.76preg20Price of energy€·kwh−10.05pxys23If 1, then farmers seek economic optimum else they seek agronomic optimumdmnl1pycm17, 24Potential marginal yield per cubic meter of waterkg·m−330sfmz19Average salinity factor across a line segment connecting the FB and ABkg·m−3impslcl8, 18, 19Crop-limiting salinitykg·m−32slfp18, 19Crop-specific calibration parameterkg·m−30.5slfw7, 8, 19Freshwater salinitykg·m−30.15slmz19Actual average salinity across a line segment connecting FB and AB, CDAB ≤x ≤ CDFBkg·m−3impslsw7, 8Saltwater salinitykg·m−335subh21Subsidies per hectare€·ha−1·year−10tfex15, 22, 25Time for farmers to adjust expectationsyear3tmbf6Time for current SB and FB to achieve equilibrium, freshening phaseyear10tmbi6Time for current SB and FB to achieve equilibrium, intrusion phaseyear5tmhe10Time for the FGS to achieve equilibriumyear5vtmz4Vertical thickness of the MZm20wdth1, 11, 26Coastal aquifer widthkm4wrao17, 23, 24, 25Water requirement for agronomic optimum yieldm3·ha−1·year−17500subscript 0Time zero, initial state of the systemyearuntsubscript sTime stepyear0.0078125subscript tState of the system at a given timeyearunt1 – dmnl = dimensionless; ha = hectare, 104 m2; € = EU Euros; 2 – imp = intermediate modelling parameter; unt = undefined time start or span.


### 3.2. Net Aquifer Discharge

One of the main hypotheses of the AQUACOAST model is that fresh groundwater flow moves towards the sea and discharges at the shoreline in natural regime. Net aquifer discharge (TDFW) is expressed as difference in inflow and outflow water balance components and is measured per unit thickness of the aquifer. In a disturbed regime, the model assumes that TDFW decreases due to groundwater pumping and the SFI displaces inland. The depth and extent of the SFI determine the fraction of usable fresh groundwater. The time-zero TDFW responds to initial pumping. This initial TDFW can be used to set other initial states in the AQUACOAST model. TDFW is expressed as(1a)TDFWt=appt⋅1−cfro−aaet−PUMPt⋅1−cfrf⋅lgth365+lttf⋅103wdth⋅365,(1b)TDFW0=appt⋅1−cfro−aaet−PUMP0⋅1−cfrf⋅lgth365+lttf⋅103wdth⋅365,where TDFW = net aquifer discharge, PUMP = groundwater pumping, appt = annual precipitation, cfro = runoff coefficient, aaet = annual actual evapotranspiration, cfrf = return flow coefficient, lgth = coastal aquifer length, lttf = lateral groundwater transference from others aquifers, wdth = coastal aquifer width, t = state of the system at a given time, and 0 = time zero.


### 3.3. Saltwater and Freshwater Boundaries Dynamics

Verruijt [39] arrives at the formulae giving the depth of the SFI below sea level (hS) and the height of the free groundwater surface (FGS: water table) (hF) as functions of the shoreline distance along x-axis, TDFW, and three constants: hydraulic conductivity (cfpm), freshwater density (dtfw), and saltwater density (dtsw). Note that x, hS, and hF, as such, are not variables of the AQUACOAST model despite that they are included in Figure 1 for illustrative purposes.As a convenient simplification, the model assumes that the depth of all the pumping wells below sea level (dpwl) is the same (Figure1). In real aquifers, pumping wells may show substantial differences in depth. In this case, dpwl will adopt the depth of the deepest one or the depth of other shallowest ones as the SFI rises due to progressive pumping.Let us call salinization boundary (SB) the line, in theory parallel to the shoreline for an isotropic and hydraulically homogeneous aquifer (whose projection in Figure1 is a point) where the SFI intersects the horizontal plane that connects the bottom of pumping wells (whose projection in Figure 1 is a line). It may be thought that the distance of the SB from the shoreline could be obtained by equaling Verruijt’s hS to dpwl and solving for x. However, this would entail assuming that the SFI moves instantaneously, which is not the case in real aquifers. Therefore, it is assumed that the expression obtained as mentioned gives the equilibrium distance from the shoreline that the SB would reach given the TDFW (EDSB) as in Figure 1. This expression includes a MIN function that avoids EDSB to be greater than the coastal aquifer length, as follows:(2a)EDSBt=MINlgth,dpwl2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFWt−TDFWt⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,(2b)EDSB0=MINlgth,dpwl2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFWt−TDFW0⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,(2c)dtcp=dtsw–dtfwdtfw,where EDSB = equilibrium distance of the SB from the shoreline, TDFW = net aquifer discharge, lgth = coastal aquifer length, dpwl = depth of pumping wells below sea level, cfpm = hydraulic conductivity, dtcp = density contrast parameter, dtsw = saltwater density, dtfw = freshwater density, t = state of the system at a given time, and 0 = time zero.The AQUACOAST model assumes that the SB moves towards equilibrium following an exponential fitting. Thus, the current distance of the SB from the shoreline (CDSB) is given by a first-order linear negative feedback. As deduced, no pumping wells whose distance from the shoreline is less than CDSB take saltwater (Figure1). It is assumed that, initially, the system is in the equilibrium associated with initial TDFW as follows:(3a)CDSBt+s=CDSBt+s⋅EDSBt−CDSBtTMBEt,(3b)CDSB0=EDSB0,where CDSB = current distance of the SB from the shoreline, EDSB = equilibrium distance of the SB from the shoreline, TMBE = time for the SB and FB to achieve equilibrium, t = state of the system at a given time, 0 = time zero, and s = time step.The AQUACOAST model represents the mixing zone (MZ) between the interface and the freshwater body (at left of the interface in Figure1). The MZ has a constant vertical thickness (vtmz) (Figure 1). Let us call freshwater boundary (FB) the line, in theory parallel to the shoreline for an isotropic and hydraulically homogeneous aquifer, where the upper-inland surface of the MZ intersects the horizontal plane connecting the bottom of pumping wells (Figure 1). The distance of the FB from the shoreline in the equilibrium associated with the current TDFW (EDFB) results from substituting dpwl by dpwl + vtmz in equations (2a)–(2c) as follows:(4a)EDFBt=MINlgth,dpwl+vtmz2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFWt−TDFWt⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,(4b)EDFB0=MINlgth,dpwl+vtmz2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFW0−TDFW0⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,where EDFB = equilibrium distance of the FB from the shoreline, TDFW = net aquifer discharge, lgth = coastal aquifer length, dpwl = depth of pumping wells below sea level, vtmz = vertical thickness of the MZ, cfpm = hydraulic conductivity, dtcp = density contrast parameter, t = state of the system at a given time, and 0 = time zero.Again, the AQUACOAST model assumes that the FB moves towards equilibrium following an exponential fitting. Thus, the current distance of the FB from the shoreline (CDFB) is given by a first-order linear negative feedback. The two surfaces limiting the MZ would move at the same velocity, so equations (3a) and (3b) for CDSB and equations (5a) and (5b) for CDFB use the same average adjustment time. As deduced, no pumping wells whose distance from the shoreline is greater than CDFB take freshwater (Figure 1). It is assumed that, initially, the system is in the equilibrium associated with the current TDFW as follows:(5a)CDFBt+s=CDFBt+s⋅EDFBt−CDFBtTMBEt,(5b)CDFB0=EDFB0,where CDFB = current distance of the FB from the shoreline, EDFB = equilibrium distance of the FB from the shoreline, TMBE = time for the SF and SB to achieve equilibrium, t = state of the system at a given time, 0 = time zero, and s = time step.The time needed for SB and FB to achieve equilibrium (TMBE) differs when saltwater is encroaching (intrusion phase) and when freshwater is rinsing the salinized aquifer (freshening phase) [40]. In real aquifers, the freshening phase is quite longer than the intrusion phase. This is the reason why saltwater encroached during the last global sea-level maximum may still remain in some coastal aquifers [4, 44]. The TMBE expression includes a conditional function to reflect this difference as(6)TMBEt=IFEDSBt>CDSBtTHENtmbi ELSE tmbf,where TMBE = time for the SB and FB to achieve equilibrium, EDSB = equilibrium distance of the SB from the shoreline, CDSB = current distance of the SB from the shoreline, tmbi = time for current boundaries to achieve equilibrium, intrusion phase, tmbf = time for current boundaries to achieve equilibrium, freshening phase, and t = state of the system at a given time.


### 3.4. Water Salinity

Verruijt’s analytical solution assumes that water salinity in the FB (at distance CDFB from the shoreline) equals freshwater salinity (slfw), whereas water salinity in the SB (at distance CDSB from the shoreline) equals saltwater salinity (slsw). The AQUACOAST model defines the water salinity as the mass fraction of total dissolved solids relative to the sample weight: pristine coastal fresh groundwater is around 0.2 kg m−3 [10, 45], fresh groundwater in yielded coastal aquifers is typically in the 0.2–0.5 kg m−3 range [46], saltwater is 35 kg m−3 [47], and crop-limiting irrigation water is 2 kg m−3 [43]. Let us assume that mixing water salinity (S) along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB) connecting SB and FB (Figure 1) is given by slsw and slfw as(7)Sx=slfw+slsw−slfw⋅CDFB−xCDFB−CDSBcdap,where S = mixing water salinity along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB), slfw = freshwater salinity, slsw = saltwater salinity, CDFB = current distance of the FB from the shoreline, CDSB = current distance of the SB from the shoreline, and cdap = calibration parameter to represent nonlinear variations in salinity.Once defined S and slfw, the distance of the abandonment boundary (AB) (CDAB) from the shoreline (Figure1) is obtained by equaling equation 7 to the crop-limiting salinity (slcl), i.e., the maximum water salinity for crop production, and then solving for x. Note that no pumping well whose distance from the shoreline is less than CDAB (i.e., slfw ˃ slcl) can use groundwater for irrigation (Figure 1). CDAB is expressed as(8)CDABt=CDFBt−CDFBt−CDSBt⋅slcl−slfwslsw−slfw1/cdap,where CDAB = current distance of the AB from the shoreline, CDFB = current distance of the SB from the shoreline, CDSB = current distance of the SB from the shoreline, slcl = crop-limiting salinity, slfw = freshwater salinity, slsw = saltwater salinity, cdap = calibration parameter to represent nonlinear variations in salinity, and t = state of the system at a given time.


### 3.5. Free Groundwater Surface

The FGS (water table) is not a horizontal plane in Verruijt’s conceptualization (Figure1). However, the AQUACOAST model only makes use of its average height above sea level. The expression giving such an average FGS in the equilibrium associated with the current TDFW (EHFW) results from integrating the water table height above sea level between 0 and the coastal aquifer length (lgth) and then dividing the result by lgth as(9a)EHFWt=8⋅dtcp⋅TDFWt⋅lgth⋅1039⋅cfpm⋅1+dtcp1/2,(9b)EHFW0=8⋅dtcp⋅TDFW0⋅lgth⋅1039⋅cfpm⋅1+dtcp1/2,where EHFW = equilibrium average height of the FGS above sea level, TDFW = net aquifer discharge, lgth = coastal aquifer length, cfpm = hydraulic conductivity, t = state of the system at a given time, and 0 = time zero.Once again, the current average height of the FGS (CHFW) is assumed to move towards equilibrium following an exponential fitting as(10a)CHFWt+s=CHFWt+sEHFWt−CHFWttmhe,(10b)CHFW0=EHFW0,where CHFW = current average height of the FGS above sea level, EHFW = equilibrium average height of the FGS above sea level, tmhe = time for the FGS to achieve equilibrium, t = state of the system at a given time, 0 = time zero, and s = time step.


### 3.6. Irrigated Area

The potential irrigated area (ARPI) is the land area over the coastal aquifer where a pumping well with depth dpwl below sea level would take fresh groundwater whose salinity does not exceed the crop-limiting salinity, i.e., the area covered by distance lgth-CDAB (Figure1). ARPI is expressed as(11)ARPIt=wdth⋅lgth−CDABt⋅102,where ARPI = potential irrigated area, CDAB = current distance of the AB from the shoreline, wdth = coastal aquifer width, lgth = coastal aquifer length, hekm = hectare-to-square-kilometer conversion parameter, and t = state of the system at a given time.Let us assume that a number of people own equally sized pieces of agricultural land over the coastal aquifer. These persons consider entering or leaving irrigated farming by comparing the average expected profit per irrigated farm (EXPF) with the returns from other alternative economic activities, i.e., with their respective opportunity costs. However, it takes some time for them to decide and carry out their plans. Additionally, suppose that the opportunity cost follows a generalized Rayleigh distribution across farmers, which is commonly used to analyze skewed data [48]. The cumulative distribution function of this distribution is generically expressed as(12)PX≤x=1−EXP−x2A2C2+1.In equation (12), X is the opportunity cost of a farmer and A and C are the mean and the coefficient of variation of the opportunity cost across farmers, respectively. Note that P, X, x, A, and C, such as, are not variables of the AQUACOAST model. In this way, P(X ≤ EXPF) is the fraction of owners whose opportunity cost is less than the average expected profit per irrigated farm, i.e., the fraction of owners that regard irrigation as profitable. But, given that farms are equally sized, P(X ≤ EXPF) is also the fraction of the potential irrigated area where irrigation is seen as profitable. Multiplying this fraction by the potential irrigated area, the target for the irrigated area (ARTI) is obtained as(13)ARTIt=ARPIt⋅1−EXP−EXPFt2ocfm2⋅ocfv2+1,where ARTI = target for the irrigated area, ARPI = potential irrigated area, EXPF = expected profit per irrigated farm, ocfm = average opportunity cost for a farmer, ocfv = coefficient of variation of the opportunity cost for a farmer, and t = state of the system at a given time.If irrigated farms are uniformly distributed over the coastal aquifer (excluding the coastal fringe CDAB where slfw ˃ slcl), the rate of abandonment of irrigated area (RAIA), because pumping wells take groundwater whose salinity exceeds the crop-limiting salinity, equals the rate of saltwater encroachment. Thus, the expression of RAIA results from dividing the velocity of the SB displacement inland by the length of the potential irrigated area. RAIA will equal zero in the freshening phase, i.e., when EDSB < CDSB. The RAIA expression is(14)RAIAt=MAX0,EDSBt−CDSBtTMBEt⋅lgth−CDABt,where RAIA = rate of abandonment of irrigated area, EDSB = equilibrium distance of the SB from the shoreline, CDSB = current distance of the SB from the shoreline, TMBE = time for the SB and FB to achieve equilibrium, CDAB = current distance of the AB from the shoreline, lgth = coastal aquifer length, and t = state of the system at a given time.The current irrigated area (ARIR) approaches towards the target for the irrigated area (ARTI) following an exponential fitting though some irrigated area is simultaneously abandoned when saltwater is encroaching, i.e., when RAIA > 0. The AQUACOAST model assumes that there is only one irrigated farm pumping from the aquifer at time zero. ARIR is expressed as(15a)ARIRt+s=ARIRt+s⋅ARTIt−ARIRttfex−RAIAt⋅ARIRt,(15b)ARIR0=arfr,where ARIR = current irrigated area, ARTI = target for the irrigated area, RAIA = rate of abandonment of irrigated area, tefx = time for farmers to adjust expectations, arfr = area of one farm, t = state of the system at a given time, and 0 = time zero.


### 3.7. Profit per Irrigated Farm

The parameter oaag allows the user to define whether the AQUACOAST model represents an open-air (oaag = 1) or an under plastic agricultural system (oaag ≠ 1). In the former case, precipitation contributes to irrigation, whereas, in the latter, it does not. The contribution of precipitation to irrigation (PPIR) is expressed as(16)PPIRt=IF oaag=1THENappt⋅1−cfro⋅10ELSE0,where PPIR = contribution of precipitation to irrigation; oaag = if 1, then open-air agriculture else under plastic agriculture; appt = annual precipitation; cfro = runoff coefficient; mmhe = cubic-meter-to-millimeter-per-hectare conversion parameter; and t = state of the system at a given time.The crop yield per hectare (YLDH) is related to the sum of the irrigation dose from groundwater (IRDS) and the contribution of precipitation to irrigation (PPIR) by means of a logistic or parabolic function. The function involves a salinity factor that measures the negative effect of water salinity on crop yield and ranges between zero and one. It can be proved that, for any nonnull value of the salinity factor, crop yield is maximized when IRDS + PPIR equals the water requirement for agronomic optimum yield (wrao) as(17)YLDHt=pycm⋅SLFTt⋅IRDSt+PPIRt·MAX0,1−IRDSt+PPIRt2⋅wrao,where YLDH = crop yield per surface hectare, SLFT = salinity factor, IRDS = irrigation dose from groundwater, PPIR = contribution of precipitation to irrigation, pycm = potential marginal yield per cubic meter of water, wrao = water requirement for agronomic optimum yield, and t = state of the system at a given time.The expression of the salinity factor (SLFT) combines two values of another salinity factor, the one given by [49]. This factor is generically expressed as(18)SLFT=MAX0,1−MAX0,S−slfpslcl−slfp,where SLFT = salinity factor generically expressed, S = mixing water salinity along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB) as in equation 7, slcl = crop-limiting salinity, slfp = another crop-specific calibration parameter, and t = state of the system at a given time.The first value of SLFT is one, its length is lgth–CDFB (Figure1), and it corresponds to the freshwater zone where pumping wells take fresh groundwater. The second value of SLFT corresponds to the inland part of the MZ where pumping wells take salty groundwater still useful for crop production (slcl ≥ slfw), and its length is CDFB–CDAB (Figure 1). Integrating equation 7 between CDAB and CDFB and then dividing the result by CDFB–CDAB gives the actual average salinity across the line segment connecting FB and AB (slmz). Thus, the value of SLFT in the second zone (sfmz) results from using slmz instead of S in equation (19). Finally, SLFT is a weighted average of the two mentioned values of SLFT, where the weights are the lengths of the respective zones. Note that, in doing so, (1) the crop yield per hectare (YLDH) is also a weighted average of the crop yields corresponding to both zones, and (2) if salinity increased linearly along the line segment between CDAB and CDFB, i.e., if cdap = 1, slmz would be the simple average of the crop limiting (slcl) and the freshwater (slfw) salinities. SLFT is expressed as(19a)SLFTt=lgth−CDFBt+sfmz⋅CDFBt−CDABtlgth−CDABt,(19b)sfmz=MAX0,1−MAX0,slmz−slfpslcl−slfp,(19c)slmz=slfw+slcl−slfwcdap+1,where YLDH = salinity factor, CDFB = current distance of the FB from the shoreline, CDAB = current distance of the AB from the shoreline, lgth = coastal aquifer length, sfmz = average salinity factor across a line segment connecting the FB and AB, slmz = actual water salinity across a line segment connecting the FB and AB (CDAB ≤ x ≤ CDFB), slfp = calibration parameter to represent nonlinear variations in salinity, slcl = crop-limiting salinity, slfw = freshwater salinity, and t = state of the system at a given time.The cost of water per cubic meter (CSTW) combines the cost of the energy needed to pump groundwater, which depends on the FGS (water table) depth and the price of energy and other costs as(20)CSTWt=egcm⋅hght−CHFWt⋅preg+cwot,where CSTW = water cost per cubic meter, CHFW = current average height of the FGS above sea level, egcm = energy required to pump one cubic meter one meter, hght = average height of land surface above sea level, preg = price of energy, cwot = costs of water other than energy, and t = state of the system at a given time.The profit per hectare (PFTH) is the difference in incomes and costs per hectare. Incomes come from the selling of agricultural production and from subsidies. The costs per hectare are the cost of water and other costs. The former results from multiplying the water cost per cubic meter (CSTW) and the actual irrigation dose from groundwater (IRDS). This product tries to reach the desired irrigation dose over the efficiency of the irrigation system. Other costs per hectare are aggregated into a single model parameter (octh). PFTH is expressed as(21)PFTHt=prcp⋅YLDHt+subh−CSTWt⋅IRDStefir−octh,where PFTH = profit per hectare, YLDH = crop yield per hectare, CSTW = water cost per cubic meter, IRDS = irrigation dose from groundwater, prcp = crop price, subh = subsidies per hectare, efir = irrigation system efficiency, octh = other costs per hectare, and t = state of the system at a given time.The expected profit per irrigated farm (EXPF) results from multiplying the profit per hectare (PFTH) and the area of one farm (arfr). The AQUACOAST model implements the expectation about such a profit is a moving average of observed past profits where weights exponentially decrease over time. This is the well-known exponential smoothing, which can be expressed as a first-order linear negative feedback [26] as(22a)EXPFt+s=EXPFt+s⋅arfr⋅PFTHt−EXPFtfex,(22b)EXPF0=expf,where EXPF = expected profit per irrigated farm, PFTH = profit per hectare, arfr = area of one farm, tfex = time for farmers to adjust expectations, expf = initial expected profit per irrigated farm, t = state of the system at a given time, 0 = time zero, and s = time step.


### 3.8. Groundwater Pumping

The AQUACOAST model includes the parameter pxys that allows the user to define whether farmers seek to maximize the profit per irrigated farm (pxys = 1) or crop yield (pxys ≠ 1). Note that the water requirement for each of these goals is different. In the latter case, the water requirement is given by the parameter wrao as in equation (17). The expression of the expected water requirement for profit maximization (EXWX) will be given later. Whatever the case, the irrigation dose from groundwater (IRDS) is the difference in water requirement and contribution of precipitation to irrigation. The resultant irrigation dose cannot be greater than the irrigation quota, if any. The expression of the initial value of IRDS simply considers EXWX = wrao. IRDS is expressed as(23a)IRDSt=MINidqt,MAX0,IF pxys=1THENEXWXtELSE wrao−PPIRt,(23b)IRDS0=MINidqt,MAX0,wrao−PPIRt,where IRDS = irrigation dose from groundwater; PPIR = contribution of precipitation to irrigation; EXWX = expected water requirement for profit maximization; idqt = irrigation quota; pxys = if 1, then farmers seek economic optimum else they seek agronomic optimum; wrao = water requirement for agronomic optimum yield; t = state of the system at a given time; and 0 = time zero.The water requirement for profit maximization (WRXP) results from calculating the partial derivative of equation (21) to deduce profit per hectare with respect to the irrigation dose, equaling it to zero, and then solving for the irrigation dose as(24)WRXPt=wrao⋅1−CSTWtefir⋅prcp⋅pycm⋅SLFTt,where WRXP = water requirement for profit maximization, CSTW = water cost per cubic meter, SLFT = salinity factor, wrao = water requirement for agronomic optimum yield, efir = irrigation system efficiency, prcp = crop price, pycm = potential marginal yield per cubic meter of water, and t = state of the system at a given time.For the current year, farmers cannot know the water requirement for profit maximization because many of the factors involved are unknown a priori, e.g., crop and energy prices, or the incidence of salinization. However, farmers can form an approximation of such a water requirement for the current year (EXWX) by learning from past experiences. This learning process would follow an exponential fitting to the true water requirement for profit maximization. In the AQUACOAST model, the initial value of EXWX is taken to be the water requirement for agronomic optimum yield, a technical parameter which is easier to grasp by farmers. EXWX is expressed as(25a)EXWXt+s=EXWXt+s⋅WRXPt−EXWXttfex,(25b)EXWX0=wrao,where EXWX = expected water requirement for profit maximization, WRXP = water requirement for profit maximization, tfex = time for farmers to adjust expectations, wrao = water requirement for agronomic optimum yield, t = state of the system at a given time, 0 = time zero, and s = time step.Groundwater pumping (PUMP) results from multiplying the current irrigated area (ARIR) and the actual irrigation dose (IRDS) and then dividing by the coastal aquifer area as(26a)PUMPt=ARIRt·IRDStwdth⋅lgth⋅efir⋅103,(26b)PUMP0=ARIR0⋅IRDS0wdth⋅lgth⋅efir⋅103,where PUMP = groundwater allotment for irrigation, ARIR = current irrigated area, IRDS = irrigation dose from groundwater, wdth = coastal aquifer width, lgth = coastal aquifer length, efir = irrigation system efficiency, t = state of the system at a given time, and 0 = time zero.


## 3.1. Overall Description

The AQUACOAST model implements the two-dimensional saltwater-freshwater interface (SFI) analytical solution proposed in [39] for an idealized unconfined (i.e., its upper boundary is the fluctuating water table) coastal aquifer, which is horizontal, enough thick, and overlying an impervious bedrock, isotropic, and hydraulically homogeneous with steady inflow (aquifer recharge) and outflow (aquifer discharge) rates (Figure 1). The variables and parameters of the AQUACOAST model concerning irrigation practices, groundwater hydrology, and groundwater allotment for agriculture have been added to the schematization of Verruijt [39]. In the steady schema of Verruijt [39], the SFI position is reached automatically depending on the fresh groundwater flow. AQUACOAST modifies this hypothesis by introducing delays in order to reach the equilibrium described by Verruijt, which magnitude differs, which differ for the intrusion phase (minor delay) and the freshening phase (longer delay) [40]. Groundwater allotment for agricultural use is the main flow component. The model focuses on determining critical distances of agronomic interest from the shoreline, such as the current and equilibrium freshwater and salinization boundaries or the irrigation farming abandonment boundary (Figure 1).Figure 1
The AQUACOAST model added to Verruijt’s schematization for an idealized unconfined coastal aquifer.hF is the height of the free groundwater surface (FGS: water table) above sea level (it is not a model parameter), hS is the depth of the saltwater-freshwater interface (SFI) below sea level (it is not a model parameter), hght is the average height of land surface above sea level, vtmz is the vertical thickness of the mixing zone, dpwl is the depth of pumping wells below sea level, CDSB is the current distance of the SB from the shoreline, CDAB is the current distance of the AB from the shoreline, CDFB is the current distance of the FB from the shoreline, EDSB is the equilibrium distance of the SB from the shoreline, EDFB is the equilibrium distance of the FB from the shoreline, lgth is the coastal aquifer length, slfw is the freshwater salinity, S is the mixing water salinity along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB), slmz is the actual average salinity across a line segment connecting FB and AB (CDAB ≤ x ≤ CDFB), slcl is the crop-limiting salinity, and slsw is the saltwater salinity.AQUACOAST is a lumped spatial model [27], and thus, its variables represent either totals or averages over the entire area modelled. For the sake of illustration, the simulation period is 50 years, an appropriate time horizon when studying environmental degradation that can only be detected when observed the evolution of “slow variables” [41]. Variables and parameters use the SI units, capital letters being variables, whereas small letters denote parameters. Tables 1 and 2 include the notation, units, and adopted values for variables and parameters, respectively, after data from [4, 8, 20, 21, 40, 42, 43] and references therein. The AQUACOAST model uses 26 straightforward formulations divided into seven sections as follows.Table 1
Notation for variables used.
VariableEquationDefinitionUnits1ARIR15, 26Current irrigated aream−2·day−1ARPI11, 13Potential irrigated aream3·ha−1·year−1ARTI13, 15Target for the irrigated area€·year−1CDAB8, 11, 14, 19Current distance of the AB from the shorelinem3·ha−1·year−1CDFB5, 7, 8, 19Current distance of the FB from the shoreline€·ha−1·year−1CDSB3, 6, 7, 8, 14Current distance of the SB from the shoreline€·m−3CHFW10, 20Current average height of the FGS above sea levelmCSTW20, 21, 24Water cost per cubic metermEDFB4, 5Equilibrium distance of the FB from the shorelinem3·ha−1·year−1EDSB2, 3, 6, 14Equilibrium distance of the SB from the shorelinem3·ha−1·year−1EHFW9, 10Equilibrium average height of the FGS above sea levelkmEXPF13, 22Expected profit per irrigated farmkmEXWX23, 25Expected water requirement for profit maximizationkmIRDS17, 21, 23, 26Irrigation dose from groundwaterkmPFTH21, 22Profit per hectarekmPPIR16, 17, 23Contribution of direct precipitation to irrigationyear−1PUMP1, 26Groundwater pumpingdmnlRAIA14, 15Rate of abandonment of irrigated areahaS(x)7, 18Salinity of mixing water along thex-axis distance from the shoreline, CDSB ≤ x ≤ CDFBkg·m−3SLFT17, 18, 19, 24Salinity factorhaTDFW1, 2, 4, 9Net aquifer dischargehaTMBE3, 5, 6, 14Time for the SB and FB to achieve equilibriumyearWRXP24, 25Water requirement for profit maximizationmm·year−1YLDH17, 21Crop yield per surface unitkg·ha−1·year−11 – dmnl = dimensionless; ha = hectare; € = EU Euros.Table 2
Notation for parameters used.
ParameterEquationDefinitionUnits1Value2aaet1Annual actual evapotranspirationmm·year−1200appt1, 16Annual precipitationmm·year−1300arfr15, 22Area of one farmha2.84cdap7, 8, 19Calibration parameter to represent nonlinear variations in salinitydmnl3cfpm2, 4, 9Hydraulic conductivitym·day−110cfrf1Return flow coefficientdmnl0.2cfro1, 16Runoff coefficientdmnl0.2cwot20Costs of water other than energy€·m−31.4dpwl2, 4Depth of pumping wells below sea levelm100dtcp2, 4, 9Density contrast parameterdmnl0.025dtfw2, 4, 9Freshwater densitykg·m−31.000dtsw2, 4, 9Saltwater densitykg·m−31.025efir21, 24, 26Irrigation system efficiencydmnl0.8egcm20Energy required to pump one cubic meter one meterkwh·m−40.017expf22Initial expected profit per irrigated farm€·year−123,000hght20Average height of land surface above sea levelm100idqt23Irrigation quotam3·ha−1·year−120,000lgth1, 2, 4, 9, 11, 14, 19, 26Coastal aquifer lengthkm10lttf1Lateral transference from other aquifershm3·year−17.9oaag16If 1, then open-air agriculture else under plastic agriculturedmnl0ocfm13Average opportunity cost for a farmer€·yr−123,000ocfv13Coefficient of variation of the opportunity cost for a farmerdmnl0.29octh21Other costs per hectare€·ha−1·year−140,053prcp21, 24Crop price€·kg−10.76preg20Price of energy€·kwh−10.05pxys23If 1, then farmers seek economic optimum else they seek agronomic optimumdmnl1pycm17, 24Potential marginal yield per cubic meter of waterkg·m−330sfmz19Average salinity factor across a line segment connecting the FB and ABkg·m−3impslcl8, 18, 19Crop-limiting salinitykg·m−32slfp18, 19Crop-specific calibration parameterkg·m−30.5slfw7, 8, 19Freshwater salinitykg·m−30.15slmz19Actual average salinity across a line segment connecting FB and AB, CDAB ≤x ≤ CDFBkg·m−3impslsw7, 8Saltwater salinitykg·m−335subh21Subsidies per hectare€·ha−1·year−10tfex15, 22, 25Time for farmers to adjust expectationsyear3tmbf6Time for current SB and FB to achieve equilibrium, freshening phaseyear10tmbi6Time for current SB and FB to achieve equilibrium, intrusion phaseyear5tmhe10Time for the FGS to achieve equilibriumyear5vtmz4Vertical thickness of the MZm20wdth1, 11, 26Coastal aquifer widthkm4wrao17, 23, 24, 25Water requirement for agronomic optimum yieldm3·ha−1·year−17500subscript 0Time zero, initial state of the systemyearuntsubscript sTime stepyear0.0078125subscript tState of the system at a given timeyearunt1 – dmnl = dimensionless; ha = hectare, 104 m2; € = EU Euros; 2 – imp = intermediate modelling parameter; unt = undefined time start or span.


## 3.2. Net Aquifer Discharge

One of the main hypotheses of the AQUACOAST model is that fresh groundwater flow moves towards the sea and discharges at the shoreline in natural regime. Net aquifer discharge (TDFW) is expressed as difference in inflow and outflow water balance components and is measured per unit thickness of the aquifer. In a disturbed regime, the model assumes that TDFW decreases due to groundwater pumping and the SFI displaces inland. The depth and extent of the SFI determine the fraction of usable fresh groundwater. The time-zero TDFW responds to initial pumping. This initial TDFW can be used to set other initial states in the AQUACOAST model. TDFW is expressed as(1a)TDFWt=appt⋅1−cfro−aaet−PUMPt⋅1−cfrf⋅lgth365+lttf⋅103wdth⋅365,(1b)TDFW0=appt⋅1−cfro−aaet−PUMP0⋅1−cfrf⋅lgth365+lttf⋅103wdth⋅365,where TDFW = net aquifer discharge, PUMP = groundwater pumping, appt = annual precipitation, cfro = runoff coefficient, aaet = annual actual evapotranspiration, cfrf = return flow coefficient, lgth = coastal aquifer length, lttf = lateral groundwater transference from others aquifers, wdth = coastal aquifer width, t = state of the system at a given time, and 0 = time zero.


## 3.3. Saltwater and Freshwater Boundaries Dynamics

Verruijt [39] arrives at the formulae giving the depth of the SFI below sea level (hS) and the height of the free groundwater surface (FGS: water table) (hF) as functions of the shoreline distance along x-axis, TDFW, and three constants: hydraulic conductivity (cfpm), freshwater density (dtfw), and saltwater density (dtsw). Note that x, hS, and hF, as such, are not variables of the AQUACOAST model despite that they are included in Figure 1 for illustrative purposes.As a convenient simplification, the model assumes that the depth of all the pumping wells below sea level (dpwl) is the same (Figure1). In real aquifers, pumping wells may show substantial differences in depth. In this case, dpwl will adopt the depth of the deepest one or the depth of other shallowest ones as the SFI rises due to progressive pumping.Let us call salinization boundary (SB) the line, in theory parallel to the shoreline for an isotropic and hydraulically homogeneous aquifer (whose projection in Figure1 is a point) where the SFI intersects the horizontal plane that connects the bottom of pumping wells (whose projection in Figure 1 is a line). It may be thought that the distance of the SB from the shoreline could be obtained by equaling Verruijt’s hS to dpwl and solving for x. However, this would entail assuming that the SFI moves instantaneously, which is not the case in real aquifers. Therefore, it is assumed that the expression obtained as mentioned gives the equilibrium distance from the shoreline that the SB would reach given the TDFW (EDSB) as in Figure 1. This expression includes a MIN function that avoids EDSB to be greater than the coastal aquifer length, as follows:(2a)EDSBt=MINlgth,dpwl2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFWt−TDFWt⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,(2b)EDSB0=MINlgth,dpwl2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFWt−TDFW0⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,(2c)dtcp=dtsw–dtfwdtfw,where EDSB = equilibrium distance of the SB from the shoreline, TDFW = net aquifer discharge, lgth = coastal aquifer length, dpwl = depth of pumping wells below sea level, cfpm = hydraulic conductivity, dtcp = density contrast parameter, dtsw = saltwater density, dtfw = freshwater density, t = state of the system at a given time, and 0 = time zero.The AQUACOAST model assumes that the SB moves towards equilibrium following an exponential fitting. Thus, the current distance of the SB from the shoreline (CDSB) is given by a first-order linear negative feedback. As deduced, no pumping wells whose distance from the shoreline is less than CDSB take saltwater (Figure1). It is assumed that, initially, the system is in the equilibrium associated with initial TDFW as follows:(3a)CDSBt+s=CDSBt+s⋅EDSBt−CDSBtTMBEt,(3b)CDSB0=EDSB0,where CDSB = current distance of the SB from the shoreline, EDSB = equilibrium distance of the SB from the shoreline, TMBE = time for the SB and FB to achieve equilibrium, t = state of the system at a given time, 0 = time zero, and s = time step.The AQUACOAST model represents the mixing zone (MZ) between the interface and the freshwater body (at left of the interface in Figure1). The MZ has a constant vertical thickness (vtmz) (Figure 1). Let us call freshwater boundary (FB) the line, in theory parallel to the shoreline for an isotropic and hydraulically homogeneous aquifer, where the upper-inland surface of the MZ intersects the horizontal plane connecting the bottom of pumping wells (Figure 1). The distance of the FB from the shoreline in the equilibrium associated with the current TDFW (EDFB) results from substituting dpwl by dpwl + vtmz in equations (2a)–(2c) as follows:(4a)EDFBt=MINlgth,dpwl+vtmz2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFWt−TDFWt⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,(4b)EDFB0=MINlgth,dpwl+vtmz2⋅cfpm⋅dtcp⋅1+dtcp2⋅TDFW0−TDFW0⋅1−dtcp2⋅dtcp⋅cfpm⋅10−3,where EDFB = equilibrium distance of the FB from the shoreline, TDFW = net aquifer discharge, lgth = coastal aquifer length, dpwl = depth of pumping wells below sea level, vtmz = vertical thickness of the MZ, cfpm = hydraulic conductivity, dtcp = density contrast parameter, t = state of the system at a given time, and 0 = time zero.Again, the AQUACOAST model assumes that the FB moves towards equilibrium following an exponential fitting. Thus, the current distance of the FB from the shoreline (CDFB) is given by a first-order linear negative feedback. The two surfaces limiting the MZ would move at the same velocity, so equations (3a) and (3b) for CDSB and equations (5a) and (5b) for CDFB use the same average adjustment time. As deduced, no pumping wells whose distance from the shoreline is greater than CDFB take freshwater (Figure 1). It is assumed that, initially, the system is in the equilibrium associated with the current TDFW as follows:(5a)CDFBt+s=CDFBt+s⋅EDFBt−CDFBtTMBEt,(5b)CDFB0=EDFB0,where CDFB = current distance of the FB from the shoreline, EDFB = equilibrium distance of the FB from the shoreline, TMBE = time for the SF and SB to achieve equilibrium, t = state of the system at a given time, 0 = time zero, and s = time step.The time needed for SB and FB to achieve equilibrium (TMBE) differs when saltwater is encroaching (intrusion phase) and when freshwater is rinsing the salinized aquifer (freshening phase) [40]. In real aquifers, the freshening phase is quite longer than the intrusion phase. This is the reason why saltwater encroached during the last global sea-level maximum may still remain in some coastal aquifers [4, 44]. The TMBE expression includes a conditional function to reflect this difference as(6)TMBEt=IFEDSBt>CDSBtTHENtmbi ELSE tmbf,where TMBE = time for the SB and FB to achieve equilibrium, EDSB = equilibrium distance of the SB from the shoreline, CDSB = current distance of the SB from the shoreline, tmbi = time for current boundaries to achieve equilibrium, intrusion phase, tmbf = time for current boundaries to achieve equilibrium, freshening phase, and t = state of the system at a given time.


## 3.4. Water Salinity

Verruijt’s analytical solution assumes that water salinity in the FB (at distance CDFB from the shoreline) equals freshwater salinity (slfw), whereas water salinity in the SB (at distance CDSB from the shoreline) equals saltwater salinity (slsw). The AQUACOAST model defines the water salinity as the mass fraction of total dissolved solids relative to the sample weight: pristine coastal fresh groundwater is around 0.2 kg m−3 [10, 45], fresh groundwater in yielded coastal aquifers is typically in the 0.2–0.5 kg m−3 range [46], saltwater is 35 kg m−3 [47], and crop-limiting irrigation water is 2 kg m−3 [43]. Let us assume that mixing water salinity (S) along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB) connecting SB and FB (Figure 1) is given by slsw and slfw as(7)Sx=slfw+slsw−slfw⋅CDFB−xCDFB−CDSBcdap,where S = mixing water salinity along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB), slfw = freshwater salinity, slsw = saltwater salinity, CDFB = current distance of the FB from the shoreline, CDSB = current distance of the SB from the shoreline, and cdap = calibration parameter to represent nonlinear variations in salinity.Once defined S and slfw, the distance of the abandonment boundary (AB) (CDAB) from the shoreline (Figure1) is obtained by equaling equation 7 to the crop-limiting salinity (slcl), i.e., the maximum water salinity for crop production, and then solving for x. Note that no pumping well whose distance from the shoreline is less than CDAB (i.e., slfw ˃ slcl) can use groundwater for irrigation (Figure 1). CDAB is expressed as(8)CDABt=CDFBt−CDFBt−CDSBt⋅slcl−slfwslsw−slfw1/cdap,where CDAB = current distance of the AB from the shoreline, CDFB = current distance of the SB from the shoreline, CDSB = current distance of the SB from the shoreline, slcl = crop-limiting salinity, slfw = freshwater salinity, slsw = saltwater salinity, cdap = calibration parameter to represent nonlinear variations in salinity, and t = state of the system at a given time.


## 3.5. Free Groundwater Surface

The FGS (water table) is not a horizontal plane in Verruijt’s conceptualization (Figure1). However, the AQUACOAST model only makes use of its average height above sea level. The expression giving such an average FGS in the equilibrium associated with the current TDFW (EHFW) results from integrating the water table height above sea level between 0 and the coastal aquifer length (lgth) and then dividing the result by lgth as(9a)EHFWt=8⋅dtcp⋅TDFWt⋅lgth⋅1039⋅cfpm⋅1+dtcp1/2,(9b)EHFW0=8⋅dtcp⋅TDFW0⋅lgth⋅1039⋅cfpm⋅1+dtcp1/2,where EHFW = equilibrium average height of the FGS above sea level, TDFW = net aquifer discharge, lgth = coastal aquifer length, cfpm = hydraulic conductivity, t = state of the system at a given time, and 0 = time zero.Once again, the current average height of the FGS (CHFW) is assumed to move towards equilibrium following an exponential fitting as(10a)CHFWt+s=CHFWt+sEHFWt−CHFWttmhe,(10b)CHFW0=EHFW0,where CHFW = current average height of the FGS above sea level, EHFW = equilibrium average height of the FGS above sea level, tmhe = time for the FGS to achieve equilibrium, t = state of the system at a given time, 0 = time zero, and s = time step.


## 3.6. Irrigated Area

The potential irrigated area (ARPI) is the land area over the coastal aquifer where a pumping well with depth dpwl below sea level would take fresh groundwater whose salinity does not exceed the crop-limiting salinity, i.e., the area covered by distance lgth-CDAB (Figure1). ARPI is expressed as(11)ARPIt=wdth⋅lgth−CDABt⋅102,where ARPI = potential irrigated area, CDAB = current distance of the AB from the shoreline, wdth = coastal aquifer width, lgth = coastal aquifer length, hekm = hectare-to-square-kilometer conversion parameter, and t = state of the system at a given time.Let us assume that a number of people own equally sized pieces of agricultural land over the coastal aquifer. These persons consider entering or leaving irrigated farming by comparing the average expected profit per irrigated farm (EXPF) with the returns from other alternative economic activities, i.e., with their respective opportunity costs. However, it takes some time for them to decide and carry out their plans. Additionally, suppose that the opportunity cost follows a generalized Rayleigh distribution across farmers, which is commonly used to analyze skewed data [48]. The cumulative distribution function of this distribution is generically expressed as(12)PX≤x=1−EXP−x2A2C2+1.In equation (12), X is the opportunity cost of a farmer and A and C are the mean and the coefficient of variation of the opportunity cost across farmers, respectively. Note that P, X, x, A, and C, such as, are not variables of the AQUACOAST model. In this way, P(X ≤ EXPF) is the fraction of owners whose opportunity cost is less than the average expected profit per irrigated farm, i.e., the fraction of owners that regard irrigation as profitable. But, given that farms are equally sized, P(X ≤ EXPF) is also the fraction of the potential irrigated area where irrigation is seen as profitable. Multiplying this fraction by the potential irrigated area, the target for the irrigated area (ARTI) is obtained as(13)ARTIt=ARPIt⋅1−EXP−EXPFt2ocfm2⋅ocfv2+1,where ARTI = target for the irrigated area, ARPI = potential irrigated area, EXPF = expected profit per irrigated farm, ocfm = average opportunity cost for a farmer, ocfv = coefficient of variation of the opportunity cost for a farmer, and t = state of the system at a given time.If irrigated farms are uniformly distributed over the coastal aquifer (excluding the coastal fringe CDAB where slfw ˃ slcl), the rate of abandonment of irrigated area (RAIA), because pumping wells take groundwater whose salinity exceeds the crop-limiting salinity, equals the rate of saltwater encroachment. Thus, the expression of RAIA results from dividing the velocity of the SB displacement inland by the length of the potential irrigated area. RAIA will equal zero in the freshening phase, i.e., when EDSB < CDSB. The RAIA expression is(14)RAIAt=MAX0,EDSBt−CDSBtTMBEt⋅lgth−CDABt,where RAIA = rate of abandonment of irrigated area, EDSB = equilibrium distance of the SB from the shoreline, CDSB = current distance of the SB from the shoreline, TMBE = time for the SB and FB to achieve equilibrium, CDAB = current distance of the AB from the shoreline, lgth = coastal aquifer length, and t = state of the system at a given time.The current irrigated area (ARIR) approaches towards the target for the irrigated area (ARTI) following an exponential fitting though some irrigated area is simultaneously abandoned when saltwater is encroaching, i.e., when RAIA > 0. The AQUACOAST model assumes that there is only one irrigated farm pumping from the aquifer at time zero. ARIR is expressed as(15a)ARIRt+s=ARIRt+s⋅ARTIt−ARIRttfex−RAIAt⋅ARIRt,(15b)ARIR0=arfr,where ARIR = current irrigated area, ARTI = target for the irrigated area, RAIA = rate of abandonment of irrigated area, tefx = time for farmers to adjust expectations, arfr = area of one farm, t = state of the system at a given time, and 0 = time zero.


## 3.7. Profit per Irrigated Farm

The parameter oaag allows the user to define whether the AQUACOAST model represents an open-air (oaag = 1) or an under plastic agricultural system (oaag ≠ 1). In the former case, precipitation contributes to irrigation, whereas, in the latter, it does not. The contribution of precipitation to irrigation (PPIR) is expressed as(16)PPIRt=IF oaag=1THENappt⋅1−cfro⋅10ELSE0,where PPIR = contribution of precipitation to irrigation; oaag = if 1, then open-air agriculture else under plastic agriculture; appt = annual precipitation; cfro = runoff coefficient; mmhe = cubic-meter-to-millimeter-per-hectare conversion parameter; and t = state of the system at a given time.The crop yield per hectare (YLDH) is related to the sum of the irrigation dose from groundwater (IRDS) and the contribution of precipitation to irrigation (PPIR) by means of a logistic or parabolic function. The function involves a salinity factor that measures the negative effect of water salinity on crop yield and ranges between zero and one. It can be proved that, for any nonnull value of the salinity factor, crop yield is maximized when IRDS + PPIR equals the water requirement for agronomic optimum yield (wrao) as(17)YLDHt=pycm⋅SLFTt⋅IRDSt+PPIRt·MAX0,1−IRDSt+PPIRt2⋅wrao,where YLDH = crop yield per surface hectare, SLFT = salinity factor, IRDS = irrigation dose from groundwater, PPIR = contribution of precipitation to irrigation, pycm = potential marginal yield per cubic meter of water, wrao = water requirement for agronomic optimum yield, and t = state of the system at a given time.The expression of the salinity factor (SLFT) combines two values of another salinity factor, the one given by [49]. This factor is generically expressed as(18)SLFT=MAX0,1−MAX0,S−slfpslcl−slfp,where SLFT = salinity factor generically expressed, S = mixing water salinity along the x-axis distance from the shoreline (CDSB ≤ x ≤ CDFB) as in equation 7, slcl = crop-limiting salinity, slfp = another crop-specific calibration parameter, and t = state of the system at a given time.The first value of SLFT is one, its length is lgth–CDFB (Figure1), and it corresponds to the freshwater zone where pumping wells take fresh groundwater. The second value of SLFT corresponds to the inland part of the MZ where pumping wells take salty groundwater still useful for crop production (slcl ≥ slfw), and its length is CDFB–CDAB (Figure 1). Integrating equation 7 between CDAB and CDFB and then dividing the result by CDFB–CDAB gives the actual average salinity across the line segment connecting FB and AB (slmz). Thus, the value of SLFT in the second zone (sfmz) results from using slmz instead of S in equation (19). Finally, SLFT is a weighted average of the two mentioned values of SLFT, where the weights are the lengths of the respective zones. Note that, in doing so, (1) the crop yield per hectare (YLDH) is also a weighted average of the crop yields corresponding to both zones, and (2) if salinity increased linearly along the line segment between CDAB and CDFB, i.e., if cdap = 1, slmz would be the simple average of the crop limiting (slcl) and the freshwater (slfw) salinities. SLFT is expressed as(19a)SLFTt=lgth−CDFBt+sfmz⋅CDFBt−CDABtlgth−CDABt,(19b)sfmz=MAX0,1−MAX0,slmz−slfpslcl−slfp,(19c)slmz=slfw+slcl−slfwcdap+1,where YLDH = salinity factor, CDFB = current distance of the FB from the shoreline, CDAB = current distance of the AB from the shoreline, lgth = coastal aquifer length, sfmz = average salinity factor across a line segment connecting the FB and AB, slmz = actual water salinity across a line segment connecting the FB and AB (CDAB ≤ x ≤ CDFB), slfp = calibration parameter to represent nonlinear variations in salinity, slcl = crop-limiting salinity, slfw = freshwater salinity, and t = state of the system at a given time.The cost of water per cubic meter (CSTW) combines the cost of the energy needed to pump groundwater, which depends on the FGS (water table) depth and the price of energy and other costs as(20)CSTWt=egcm⋅hght−CHFWt⋅preg+cwot,where CSTW = water cost per cubic meter, CHFW = current average height of the FGS above sea level, egcm = energy required to pump one cubic meter one meter, hght = average height of land surface above sea level, preg = price of energy, cwot = costs of water other than energy, and t = state of the system at a given time.The profit per hectare (PFTH) is the difference in incomes and costs per hectare. Incomes come from the selling of agricultural production and from subsidies. The costs per hectare are the cost of water and other costs. The former results from multiplying the water cost per cubic meter (CSTW) and the actual irrigation dose from groundwater (IRDS). This product tries to reach the desired irrigation dose over the efficiency of the irrigation system. Other costs per hectare are aggregated into a single model parameter (octh). PFTH is expressed as(21)PFTHt=prcp⋅YLDHt+subh−CSTWt⋅IRDStefir−octh,where PFTH = profit per hectare, YLDH = crop yield per hectare, CSTW = water cost per cubic meter, IRDS = irrigation dose from groundwater, prcp = crop price, subh = subsidies per hectare, efir = irrigation system efficiency, octh = other costs per hectare, and t = state of the system at a given time.The expected profit per irrigated farm (EXPF) results from multiplying the profit per hectare (PFTH) and the area of one farm (arfr). The AQUACOAST model implements the expectation about such a profit is a moving average of observed past profits where weights exponentially decrease over time. This is the well-known exponential smoothing, which can be expressed as a first-order linear negative feedback [26] as(22a)EXPFt+s=EXPFt+s⋅arfr⋅PFTHt−EXPFtfex,(22b)EXPF0=expf,where EXPF = expected profit per irrigated farm, PFTH = profit per hectare, arfr = area of one farm, tfex = time for farmers to adjust expectations, expf = initial expected profit per irrigated farm, t = state of the system at a given time, 0 = time zero, and s = time step.


## 3.8. Groundwater Pumping

The AQUACOAST model includes the parameter pxys that allows the user to define whether farmers seek to maximize the profit per irrigated farm (pxys = 1) or crop yield (pxys ≠ 1). Note that the water requirement for each of these goals is different. In the latter case, the water requirement is given by the parameter wrao as in equation (17). The expression of the expected water requirement for profit maximization (EXWX) will be given later. Whatever the case, the irrigation dose from groundwater (IRDS) is the difference in water requirement and contribution of precipitation to irrigation. The resultant irrigation dose cannot be greater than the irrigation quota, if any. The expression of the initial value of IRDS simply considers EXWX = wrao. IRDS is expressed as(23a)IRDSt=MINidqt,MAX0,IF pxys=1THENEXWXtELSE wrao−PPIRt,(23b)IRDS0=MINidqt,MAX0,wrao−PPIRt,where IRDS = irrigation dose from groundwater; PPIR = contribution of precipitation to irrigation; EXWX = expected water requirement for profit maximization; idqt = irrigation quota; pxys = if 1, then farmers seek economic optimum else they seek agronomic optimum; wrao = water requirement for agronomic optimum yield; t = state of the system at a given time; and 0 = time zero.The water requirement for profit maximization (WRXP) results from calculating the partial derivative of equation (21) to deduce profit per hectare with respect to the irrigation dose, equaling it to zero, and then solving for the irrigation dose as(24)WRXPt=wrao⋅1−CSTWtefir⋅prcp⋅pycm⋅SLFTt,where WRXP = water requirement for profit maximization, CSTW = water cost per cubic meter, SLFT = salinity factor, wrao = water requirement for agronomic optimum yield, efir = irrigation system efficiency, prcp = crop price, pycm = potential marginal yield per cubic meter of water, and t = state of the system at a given time.For the current year, farmers cannot know the water requirement for profit maximization because many of the factors involved are unknown a priori, e.g., crop and energy prices, or the incidence of salinization. However, farmers can form an approximation of such a water requirement for the current year (EXWX) by learning from past experiences. This learning process would follow an exponential fitting to the true water requirement for profit maximization. In the AQUACOAST model, the initial value of EXWX is taken to be the water requirement for agronomic optimum yield, a technical parameter which is easier to grasp by farmers. EXWX is expressed as(25a)EXWXt+s=EXWXt+s⋅WRXPt−EXWXttfex,(25b)EXWX0=wrao,where EXWX = expected water requirement for profit maximization, WRXP = water requirement for profit maximization, tfex = time for farmers to adjust expectations, wrao = water requirement for agronomic optimum yield, t = state of the system at a given time, 0 = time zero, and s = time step.Groundwater pumping (PUMP) results from multiplying the current irrigated area (ARIR) and the actual irrigation dose (IRDS) and then dividing by the coastal aquifer area as(26a)PUMPt=ARIRt·IRDStwdth⋅lgth⋅efir⋅103,(26b)PUMP0=ARIR0⋅IRDS0wdth⋅lgth⋅efir⋅103,where PUMP = groundwater allotment for irrigation, ARIR = current irrigated area, IRDS = irrigation dose from groundwater, wdth = coastal aquifer width, lgth = coastal aquifer length, efir = irrigation system efficiency, t = state of the system at a given time, and 0 = time zero.


## 4. AQUACOAST Venapp

### 4.1. Overall Description

The Vensim® DSS extension Venapps™ was utilized to program the user-friendly interface of the above described formulations. To use AQUACOAST, it is necessary to install Vensim Reader and the folder with all the necessary files (Table3) on the computer.Table 3
Files compiled in the AQUACOAST Venapp.
FileTypeNumberDescriptionAQUACOAST.vmfSimulation model1The usual extension of a Vensim model is ∗.mdl. However, for Venapp programming, it must be saved in binary format (∗.vmf).AQUACOAST.vcdVenapp1Vensim application.Intro.vmfSimulation model1Simulation model to explore Vensim options.Intro.vcdVenapp1Vensim application.Intro.vgdVensim graph dataCustomized graphs that allow displaying more than one variable, with more than one scale, on a single graph.∗.vdfVensim data filenContains results from a simulation, including scenario values; the name of the exercise (∗) is provided by the user.Default.vdfVensim data file1Contains results from a default simulation, including scenario values.After opening the AQUACOAST file, the welcome screen will appear (Figure2(a)). Here is brief information about the model and its authors. Clicking on the “Continue” button takes the user to the main menu (Figure 2(b)). From this point, it is possible to (i) review the model structure, (ii) access the screen from which the simulation scenarios configured the model is executed, and (iii) analyze the results using different tools of Vensim. The bottom part of the screen gives the option to explore the possibilities offered by Vensim (button “Find Out More about Vensim”), which is supported by the Intro.vmf model file (Table 3) and is part of the Vensim examples, or leave the simulator (button “Exit Model Reader”).Figure 2
(a) Welcome screen of the AQUACOAST Venapp; (b) main menu screen.
(a)(b)


### 4.2. Exploring the Structure of the Model

One of the main purposes of SD is to understand how a system works. To do this, Vensim offers multiple tools. The first allows studying the model structure. For this to be practical, the model must be split into different screens or “Views” by clicking the button “Review Model Structure” as in Figure2(b) so that the user can gradually assimilate the relationship between the variables involved. Figure 3 shows one of these screens as an example. The button captions themselves explain the available utilities.Figure 3
“Coastal aquifer dynamics” view of the AQUACOAST model.[[parms resize(1),pos(50,50),size(200,200),bgcol(156)]]


### 4.3. Running the Model

The button “Simulate the Model” as in Figure2(b) enables the model to be executed. In this part of the Venapp (Figure 4), the simulation scenarios are set. For the former, the “New scenario…” button opens a dialog box to select an existing simulation or create a new one. After that, the “Make changes…” button will show all the model parameters and their default values; by modifying them, alternative scenarios will be created.Figure 4
Simulation setup screen, including setting up and simulation scenarios.These changes are not permanent, for it would be necessary to access the model through a commercial license of Vensim® DSS and change them. Therefore, if the model is simulated without any modification, the user will obtain the results contained in the Default.vdf simulation file (Table3). After setting the scenario, the model is ready to be simulated. After clicking on the button to execute Vensim “Simulate the Model,” a blank screen will appear on which the user has to click to move on to the next stage, which consists of visualizing and analyzing the results.


### 4.4. Analyzing Results

After clicking button “Analyze Simulation Results” as in Figure2(b), this stage offers three blocks of possibilities: analysis control, results, and causal tracing (Figure 5(a)). The “Analysis Control” block allows the user to manage simulations, select variables, and compare scenarios. This is an useful option when starting to accumulate simulations. In that case, it is hard to remember what changes from one scenario to another. Using the “Load, unload…” button, this option details which parameters change (and how much) between the two previously selected simulations.Figure 5
(a) “Analysis of the results” screen including three blocks: analysis control, results, and causal tracing; (b) for the “result” block, menu of available graphs; (c) temporal trends of the main variables for the default optional scenarios, in this example, “irrigated area.”The “Results” block displays predesigned graphs showing the temporal evolution of the most important variables of the model (Figure5(a)). Clicking on the button “Display a predefined...,” a list will appear with the available graphics (Figure 5(b)). In all of them, the results of the default scenario will be included, as well as the results of the simulation that appear in the first position of the list of loaded simulations, which are managed from the “Load, unload, or reorder…” button, as in Figure 5(c).The “Causal Tracing” block (Figure5(a)) enables the use of Vensim’s own analysis tools, as shown in Figures 6–8. The tree diagram tool creates a graphical representation of the structure of a model associated with the selected workbench variable; the currently active variable is displayed at the top left of the screen. The tree diagram tool provides an easy-to-read, graphical representation of the causes and uses of variables. This is a convenient way for displaying structural relationships and therefore information about a complex structure in a relatively small space.Figure 6
“Trace underlying causes with Trees” screen. Example of the underlying causes of “Crop yield per hectare” variable using trees.Figure 7
“Trace underlying causes with Graphs” screen. Example of the underlying causes of “Crop yield per hectare” variable using graphs.Figure 8
“Trace the Uses of a variable” screen. Example of the underlying causes of “Crop yield per hectare” variable.The causal structure of the model can be explored using trees or charts, whereas the uses of each variable can only be analyzed using trees by displaying the tree of variables affected by the selected one. Figure6 shows the screen that appears after clicking the button “Trace underlying causes using Trees.” At the top is the causal tree of the currently active variable. It is possible to select other variables using the “Select a new variable...” button or, interactively, by clicking on any variable in the same causal tree. In addition, the user can (i) access to the model equations through the button “Definition...” and (ii) see its temporal trajectory according to the last simulated or loaded scenario. The rest of the buttons allow the user to scroll to other screens, as their own caption indicates. In this way, all the options of using the tree diagram tool are interconnected.Clicking the button “Trace underlying causes using Graphs” (Figure5(a)), the graphical option to study the causal tree of a variable can be used (Figure 7). In this case, only first-order variables, i.e., those directly affected by the selected variable, are shown. On the right side of the screen, the temporal trajectories of the variables included in the performance equation can be observed, as well as the value of the parameters involved. The definition of the variable is also presented in the box on the left. The rest of the options allow the user to navigate between the screens related to the diagram tool.The third analysis option offered by this Venapp (Figure5(a)) is to “Trace the Uses of a variable” (Figure 8). Again, it is possible to see its definition and scroll to other analysis screens. Finally, it is important to highlight the interactivity of the tree diagram tool. Any of the outputs it generates allows the user to select a new variable by clicking on it. This option is designed to make it easier to explore and study the model, one of Vensim’s hallmarks.


## 4.1. Overall Description

The Vensim® DSS extension Venapps™ was utilized to program the user-friendly interface of the above described formulations. To use AQUACOAST, it is necessary to install Vensim Reader and the folder with all the necessary files (Table3) on the computer.Table 3
Files compiled in the AQUACOAST Venapp.
FileTypeNumberDescriptionAQUACOAST.vmfSimulation model1The usual extension of a Vensim model is ∗.mdl. However, for Venapp programming, it must be saved in binary format (∗.vmf).AQUACOAST.vcdVenapp1Vensim application.Intro.vmfSimulation model1Simulation model to explore Vensim options.Intro.vcdVenapp1Vensim application.Intro.vgdVensim graph dataCustomized graphs that allow displaying more than one variable, with more than one scale, on a single graph.∗.vdfVensim data filenContains results from a simulation, including scenario values; the name of the exercise (∗) is provided by the user.Default.vdfVensim data file1Contains results from a default simulation, including scenario values.After opening the AQUACOAST file, the welcome screen will appear (Figure2(a)). Here is brief information about the model and its authors. Clicking on the “Continue” button takes the user to the main menu (Figure 2(b)). From this point, it is possible to (i) review the model structure, (ii) access the screen from which the simulation scenarios configured the model is executed, and (iii) analyze the results using different tools of Vensim. The bottom part of the screen gives the option to explore the possibilities offered by Vensim (button “Find Out More about Vensim”), which is supported by the Intro.vmf model file (Table 3) and is part of the Vensim examples, or leave the simulator (button “Exit Model Reader”).Figure 2
(a) Welcome screen of the AQUACOAST Venapp; (b) main menu screen.
(a)(b)


## 4.2. Exploring the Structure of the Model

One of the main purposes of SD is to understand how a system works. To do this, Vensim offers multiple tools. The first allows studying the model structure. For this to be practical, the model must be split into different screens or “Views” by clicking the button “Review Model Structure” as in Figure2(b) so that the user can gradually assimilate the relationship between the variables involved. Figure 3 shows one of these screens as an example. The button captions themselves explain the available utilities.Figure 3
“Coastal aquifer dynamics” view of the AQUACOAST model.[[parms resize(1),pos(50,50),size(200,200),bgcol(156)]]


## 4.3. Running the Model

The button “Simulate the Model” as in Figure2(b) enables the model to be executed. In this part of the Venapp (Figure 4), the simulation scenarios are set. For the former, the “New scenario…” button opens a dialog box to select an existing simulation or create a new one. After that, the “Make changes…” button will show all the model parameters and their default values; by modifying them, alternative scenarios will be created.Figure 4
Simulation setup screen, including setting up and simulation scenarios.These changes are not permanent, for it would be necessary to access the model through a commercial license of Vensim® DSS and change them. Therefore, if the model is simulated without any modification, the user will obtain the results contained in the Default.vdf simulation file (Table3). After setting the scenario, the model is ready to be simulated. After clicking on the button to execute Vensim “Simulate the Model,” a blank screen will appear on which the user has to click to move on to the next stage, which consists of visualizing and analyzing the results.


## 4.4. Analyzing Results

After clicking button “Analyze Simulation Results” as in Figure2(b), this stage offers three blocks of possibilities: analysis control, results, and causal tracing (Figure 5(a)). The “Analysis Control” block allows the user to manage simulations, select variables, and compare scenarios. This is an useful option when starting to accumulate simulations. In that case, it is hard to remember what changes from one scenario to another. Using the “Load, unload…” button, this option details which parameters change (and how much) between the two previously selected simulations.Figure 5
(a) “Analysis of the results” screen including three blocks: analysis control, results, and causal tracing; (b) for the “result” block, menu of available graphs; (c) temporal trends of the main variables for the default optional scenarios, in this example, “irrigated area.”The “Results” block displays predesigned graphs showing the temporal evolution of the most important variables of the model (Figure5(a)). Clicking on the button “Display a predefined...,” a list will appear with the available graphics (Figure 5(b)). In all of them, the results of the default scenario will be included, as well as the results of the simulation that appear in the first position of the list of loaded simulations, which are managed from the “Load, unload, or reorder…” button, as in Figure 5(c).The “Causal Tracing” block (Figure5(a)) enables the use of Vensim’s own analysis tools, as shown in Figures 6–8. The tree diagram tool creates a graphical representation of the structure of a model associated with the selected workbench variable; the currently active variable is displayed at the top left of the screen. The tree diagram tool provides an easy-to-read, graphical representation of the causes and uses of variables. This is a convenient way for displaying structural relationships and therefore information about a complex structure in a relatively small space.Figure 6
“Trace underlying causes with Trees” screen. Example of the underlying causes of “Crop yield per hectare” variable using trees.Figure 7
“Trace underlying causes with Graphs” screen. Example of the underlying causes of “Crop yield per hectare” variable using graphs.Figure 8
“Trace the Uses of a variable” screen. Example of the underlying causes of “Crop yield per hectare” variable.The causal structure of the model can be explored using trees or charts, whereas the uses of each variable can only be analyzed using trees by displaying the tree of variables affected by the selected one. Figure6 shows the screen that appears after clicking the button “Trace underlying causes using Trees.” At the top is the causal tree of the currently active variable. It is possible to select other variables using the “Select a new variable...” button or, interactively, by clicking on any variable in the same causal tree. In addition, the user can (i) access to the model equations through the button “Definition...” and (ii) see its temporal trajectory according to the last simulated or loaded scenario. The rest of the buttons allow the user to scroll to other screens, as their own caption indicates. In this way, all the options of using the tree diagram tool are interconnected.Clicking the button “Trace underlying causes using Graphs” (Figure5(a)), the graphical option to study the causal tree of a variable can be used (Figure 7). In this case, only first-order variables, i.e., those directly affected by the selected variable, are shown. On the right side of the screen, the temporal trajectories of the variables included in the performance equation can be observed, as well as the value of the parameters involved. The definition of the variable is also presented in the box on the left. The rest of the options allow the user to navigate between the screens related to the diagram tool.The third analysis option offered by this Venapp (Figure5(a)) is to “Trace the Uses of a variable” (Figure 8). Again, it is possible to see its definition and scroll to other analysis screens. Finally, it is important to highlight the interactivity of the tree diagram tool. Any of the outputs it generates allows the user to select a new variable by clicking on it. This option is designed to make it easier to explore and study the model, one of Vensim’s hallmarks.


## 5. Discussion and Conclusions

### 5.1. Global Relevance of Coastal Groundwater Resources

Groundwater is one of the more strategic renewable resources, especially in coastal drylands [4]. Over two billion people rely on groundwater as their primary potable freshwater source [50], whereas half or more of the irrigation water used for food production has also this source [17]. Groundwater also acts as the key strategic reserve in times of drought [51] and to set population in drylands with sparse economic options except agriculture [8, 20].Soaring water withdrawal worsens water scarcity already prevalent in those coastal drylands having limited surface water resources. A consequence of the aquifer overexploitation to cope with water requirements is saltwater intrusion, which threats the sustainable food production and economic development in an increasing number of countries [15, 52, 53].The use of SD modelling in water management has a long tradition due to the complex nature of the problems addressed [54]. However, few models consider the interaction between terrestrial water fluxes and human activities and associated uses, and even fewer are the models intended to coastal groundwater-dependent irrigation agriculture [15]. The AQUACOAST tool is a contribution to fill this gap of knowledge.The AQUACOAST is a global application tool despite that the Mediterranean region is postulated as a main scenario in which to apply. Coastal aquifers in this area support strong groundwater pumping rates for crop production, and saltwater intrusion has extended in most of them (Figure9). For instance, in Spain, this problem affects most of the Mediterranean coastal aquifers, so they were considered one of the “desertification landscapes” in the National Program to Combat Desertification [56, 57]. Our immediate purpose is to apply AQUACOAST to coastal aquifers of the Mediterranean arc where tourist and agricultural pressure has triggered problems of marine intrusion such as Campo de Dalías, Campo de Níjar, Vélez River and Verde River in Spain, Martíl-Alila and Smir in Morocco, and Nador and Tipaza in Argelia.Figure 9
Southern Mediterranean European groundwater bodies affected by saltwater intrusion due to overexploitation (modified from [55]).


### 5.2. Using SD for Modelling Coastal Groundwater-Dependent Irrigation Systems

SD models express in mathematical language (the most precise way that exists) the dynamic processes that take place in a given field of study. The academic and scientific value of understanding how many variables relate, determining what is cause and what is effect, making explicit feedback loops of a system, and, finally, translating them into a set of operative differential equations is very high.However, it is important to make an additional effort so that this information reaches the less specialized public involved in decision-making processes. In order to do this, the use of DSSs is essential. An example taken from AQUACOAST can help us understand the value of programming an application aimed at a wider audience. Figure10 shows the causal diagram drawn up for purely academic use as an indecipherable amalgam in the eyes of someone other than the modeler. Variables are represented by acronyms, and it is difficult to follow the path of causality in the complex web of causal relationships. However, a careful edition of the model, expliciting the full name of each variable and splitting the diagram into five screens (Figure 3 shows one of them), allows the user to create a much more user-friendly version.Figure 10
Original sketch of the AQUACOAST model. Acronyms for variables (capital letters) and parameters (small letters) are as in Tables1 and 2, respectively.The precision of mathematical language should not confuse us about the nature of the SD models. In reality, when confronted with socioecological systems, we lack universal physical or chemical laws that allow us to talk about predictions. Socioecological systems require input from a wide range of sources and types of knowledge [58]. This includes qualitative and quantitative pieces on information from various stakeholder groups such as scientists, policy makers, and community members. To collect, synthesize, and use these datasets in useful ways, integrated assessment and modelling needs to utilize and combine different methods (i.e., conceptual, numerical, and participatory) in appropriate methodological designs that best fit the project’s context, objective, and constraints, in the latter case including resource availability [27, 33].The AQUACOAST model is not intended for prediction or forecasting, even though they provide outputs over variable time periods. This is because there are not enough data to validate the models for such purposes; remember, for example, the lack of some essential information such as the initial groundwater storage or similar others like past weather data or former land uses. Thus, the aim is to get qualitative rather than quantitative outputs, answering the “what if” questions, which helps to compare simulation scenarios that address different hypotheses.System understanding and experimentation also constitute important goals of the models, as well as social learning, though this is limited to those cases where the assessment shows a special risk of degradation. Since degradation usually proceeds slowly, stakeholders are not sufficiently aware of its consequences, no matter how they value the water and land use being degraded [59]. Thus, the AQUACOAST model is intended as a “means of exploration” [60, 61] for helping them to better understand how a coastal groundwater-dependent agricultural system may behave.


### 5.3. The AQUACOAST Venapp

Just as there are several SD modelling platforms, DSS can be built using a wide range of possibilities. Venapps offer clear advantages when the software chosen is VensimDSS. The construction is very fast since there are templates that with few modifications allow the user to implement the main options of a DSS: configuration of scenarios, simulation, and visualization of the results and their analysis. The use of the tool is simple since it only requires to install the free version of VensimDSS (Vensim Reader) and to copy the Vensim files provided by the modeler.The main problem is the graphical limitation of Vensim, which prevents the use of images or figures or the implementation of options such as global sensitivity analysis or risk estimation. More complex DSS based on Vensim can be built using various programming languages (e.g., Visual Basic) and platforms (e.g., Excel), such as the application designed by authors to study desertification risk in Spain, which combined some hotspots such as soil erosion, grassland overgrazing, and aquifer overexploitation [57]. However, the choice of the type of DSS depends, as does the construction of a model, on its purpose.Since the AQUACOAST tool is aimed at (i) facilitating the study and understanding of coastal groundwater-dependent irrigated agriculture systems and (ii) creating an application whose installation and use are very simple, Venapp is a perfect ally since it makes available to the user, in a very simple way, all the analysis tools of Vensim.Since Meadows et al. [62] presented the revolutionary best-seller “Limits to growth,” for which they used systems thinking and SD concepts in order to explain how short-term development policies can lead to “overshoot and collapse” behavior of socioecological systems, SD has shown its potential as a tool to help understand complex socioecological systems and is still regarded as a valuable resource for thinking about sustainable futures [33, 63]. The AQUACOAST tool is in line with these principles and aims to highlight the fragility of such a strategic renewable resource as groundwater. As long as the model describes reality with a certain accuracy, the modelling process and its outcomes can be used to improve our understanding of the problem of excessive pumping under the increasing food production needs a necessary step towards affecting sustainable and effective change [64].


## 5.1. Global Relevance of Coastal Groundwater Resources

Groundwater is one of the more strategic renewable resources, especially in coastal drylands [4]. Over two billion people rely on groundwater as their primary potable freshwater source [50], whereas half or more of the irrigation water used for food production has also this source [17]. Groundwater also acts as the key strategic reserve in times of drought [51] and to set population in drylands with sparse economic options except agriculture [8, 20].Soaring water withdrawal worsens water scarcity already prevalent in those coastal drylands having limited surface water resources. A consequence of the aquifer overexploitation to cope with water requirements is saltwater intrusion, which threats the sustainable food production and economic development in an increasing number of countries [15, 52, 53].The use of SD modelling in water management has a long tradition due to the complex nature of the problems addressed [54]. However, few models consider the interaction between terrestrial water fluxes and human activities and associated uses, and even fewer are the models intended to coastal groundwater-dependent irrigation agriculture [15]. The AQUACOAST tool is a contribution to fill this gap of knowledge.The AQUACOAST is a global application tool despite that the Mediterranean region is postulated as a main scenario in which to apply. Coastal aquifers in this area support strong groundwater pumping rates for crop production, and saltwater intrusion has extended in most of them (Figure9). For instance, in Spain, this problem affects most of the Mediterranean coastal aquifers, so they were considered one of the “desertification landscapes” in the National Program to Combat Desertification [56, 57]. Our immediate purpose is to apply AQUACOAST to coastal aquifers of the Mediterranean arc where tourist and agricultural pressure has triggered problems of marine intrusion such as Campo de Dalías, Campo de Níjar, Vélez River and Verde River in Spain, Martíl-Alila and Smir in Morocco, and Nador and Tipaza in Argelia.Figure 9
Southern Mediterranean European groundwater bodies affected by saltwater intrusion due to overexploitation (modified from [55]).


## 5.2. Using SD for Modelling Coastal Groundwater-Dependent Irrigation Systems

SD models express in mathematical language (the most precise way that exists) the dynamic processes that take place in a given field of study. The academic and scientific value of understanding how many variables relate, determining what is cause and what is effect, making explicit feedback loops of a system, and, finally, translating them into a set of operative differential equations is very high.However, it is important to make an additional effort so that this information reaches the less specialized public involved in decision-making processes. In order to do this, the use of DSSs is essential. An example taken from AQUACOAST can help us understand the value of programming an application aimed at a wider audience. Figure10 shows the causal diagram drawn up for purely academic use as an indecipherable amalgam in the eyes of someone other than the modeler. Variables are represented by acronyms, and it is difficult to follow the path of causality in the complex web of causal relationships. However, a careful edition of the model, expliciting the full name of each variable and splitting the diagram into five screens (Figure 3 shows one of them), allows the user to create a much more user-friendly version.Figure 10
Original sketch of the AQUACOAST model. Acronyms for variables (capital letters) and parameters (small letters) are as in Tables1 and 2, respectively.The precision of mathematical language should not confuse us about the nature of the SD models. In reality, when confronted with socioecological systems, we lack universal physical or chemical laws that allow us to talk about predictions. Socioecological systems require input from a wide range of sources and types of knowledge [58]. This includes qualitative and quantitative pieces on information from various stakeholder groups such as scientists, policy makers, and community members. To collect, synthesize, and use these datasets in useful ways, integrated assessment and modelling needs to utilize and combine different methods (i.e., conceptual, numerical, and participatory) in appropriate methodological designs that best fit the project’s context, objective, and constraints, in the latter case including resource availability [27, 33].The AQUACOAST model is not intended for prediction or forecasting, even though they provide outputs over variable time periods. This is because there are not enough data to validate the models for such purposes; remember, for example, the lack of some essential information such as the initial groundwater storage or similar others like past weather data or former land uses. Thus, the aim is to get qualitative rather than quantitative outputs, answering the “what if” questions, which helps to compare simulation scenarios that address different hypotheses.System understanding and experimentation also constitute important goals of the models, as well as social learning, though this is limited to those cases where the assessment shows a special risk of degradation. Since degradation usually proceeds slowly, stakeholders are not sufficiently aware of its consequences, no matter how they value the water and land use being degraded [59]. Thus, the AQUACOAST model is intended as a “means of exploration” [60, 61] for helping them to better understand how a coastal groundwater-dependent agricultural system may behave.


## 5.3. The AQUACOAST Venapp

Just as there are several SD modelling platforms, DSS can be built using a wide range of possibilities. Venapps offer clear advantages when the software chosen is VensimDSS. The construction is very fast since there are templates that with few modifications allow the user to implement the main options of a DSS: configuration of scenarios, simulation, and visualization of the results and their analysis. The use of the tool is simple since it only requires to install the free version of VensimDSS (Vensim Reader) and to copy the Vensim files provided by the modeler.The main problem is the graphical limitation of Vensim, which prevents the use of images or figures or the implementation of options such as global sensitivity analysis or risk estimation. More complex DSS based on Vensim can be built using various programming languages (e.g., Visual Basic) and platforms (e.g., Excel), such as the application designed by authors to study desertification risk in Spain, which combined some hotspots such as soil erosion, grassland overgrazing, and aquifer overexploitation [57]. However, the choice of the type of DSS depends, as does the construction of a model, on its purpose.Since the AQUACOAST tool is aimed at (i) facilitating the study and understanding of coastal groundwater-dependent irrigated agriculture systems and (ii) creating an application whose installation and use are very simple, Venapp is a perfect ally since it makes available to the user, in a very simple way, all the analysis tools of Vensim.Since Meadows et al. [62] presented the revolutionary best-seller “Limits to growth,” for which they used systems thinking and SD concepts in order to explain how short-term development policies can lead to “overshoot and collapse” behavior of socioecological systems, SD has shown its potential as a tool to help understand complex socioecological systems and is still regarded as a valuable resource for thinking about sustainable futures [33, 63]. The AQUACOAST tool is in line with these principles and aims to highlight the fragility of such a strategic renewable resource as groundwater. As long as the model describes reality with a certain accuracy, the modelling process and its outcomes can be used to improve our understanding of the problem of excessive pumping under the increasing food production needs a necessary step towards affecting sustainable and effective change [64].



---

*Source: 9092829-2020-05-14.xml*