A review of conjunctive GW-SW management by simulation–optimization tools

The conjunctive use of groundwater and surface water (GW-SW) resources is an effective approach to manage water resources systems (Payne 2005; Alimohammadi & Hosseinzadeh 2010; Pathare & Jagtap 2016; American Society of Civil Engineers (ASCE) 2017). Conjunctive use is practiced to allocate surface and groundwater to users based on qualitative and quantitative criteria while taking into account system constraints such as economic and social factors (Safavi et al. 2016; Portoghesea et al. 2021). Alluvial aquifers are a significant source of water supply in numerous hydrogeological regions extending along rivers (Rosenshein 1988; Larkin & Sharp 1992; Barlow et al. 2003). Alluvial aquifers occur adjacent to rivers and in buried channels. These aquifers have variable thickness and contain layers of sand and gravel deposited by networks of streams over geologic time (Heath 1984). Alluvial groundwater has a strong hydraulic connection with natural surface streams (rivers), surface runoff, lakes, and reservoirs. Therefore, groundwater pumping in these aquifers leads to streamflow depletion and reduction of lake and reservoir storage. Correct planning of groundwater withdrawal in alluvial aquifers requires knowledge of GW-SW interactions (Sophocleous 2002; Rossetto et al. 2019). Spatial and temporal interactions of GW-SW resources can be assessed by the integrated, dynamic coupling of GW-SW resources employing mathematical models (Ghordoyee Milan et al. 2018).

Optimal operational rules for GW-SW management have a high priority in alluvial basins, which commonly feature productive aquifers (Brookfield & Gnau 2016). System analysis techniques for water resource planning and management have achieved remarkable advances over the last few decades (Ahmadianfar et al. 2019a). System analysis techniques may be grouped into simulation, optimization, and simulation–optimization methods (Loucks et al. 1981). Coupling planning/operation models with hydrologic models (e.g., MODSIM-MODFLOW) are currently introduced as a powerful tool in integrated water resource management (Morway et al. 2016). Various optimization tools have been widely applied to tackle real-world problems in water resource science (Ahmadianfar et al. 2019b), in which the result of some novel optimizers such as the gradient-based method has been more promising than evolutionary optimization techniques such as the genetic algorithm (GA) and the differential evolution (DE) algorithm (Ahmadianfar et al. 2020). Thorough recognition of the GW-SW models and available coupling optimization tools allow the experts to choose the best combination to either solve or predict the problems in watersheds. Simulation and optimization techniques are combined to find possible alternatives for the conjunctive use of GW-SW resources (Karamouz et al. 2003).

Despite multiple publications on the application of optimization techniques in different fields of water resources such as irrigated agriculture (Bozorg-Haddad et al. 2009), reservoir operation (Shokri et al. 2013; Akbari-Alashti et al. 2014; Asgari et al. 2016), and water distribution networks (Soltanjalili et al. 2011), this evaluates simulation and optimization methods and their many applications. It is noteworthy that a few reviews regarding GW-SW models have been published, which focused on the use of partial differential equations to model the components of coupled systems (Paniconi & Putti 2015), discussed the physical and numerical alternatives of several models (Furman 2008), and covered the conceptual approaches to GW-SW interaction (Sophocleous 2002). However, the combination of GW-SW simulation models and related optimization techniques have not been comprehensively evaluated in review articles. Therefore, this work attempts to fill the following gaps:

• determine operational rules for modeling the conjunctive use of GW-SW resources,

• present an outlook of the current challenges and trends in the conjunctive use of GW-SW resources, and

• overview optimization methods and their applications in coupled simulation–optimization of GW-SW resources management.

This work reviews conjunctive GW-SW models coupled with optimizers; yet the large volume of existing publications does not allow a complete coverage of all published works, but, rather, a selective review of the main types of approaches following in the subject matter of GW-SW simulation–optimization.

A systematic understanding of GW-SW interactions is necessary to manage scarce water resources in alluvial aquifers in arid and semi-arid regions (Wu et al. 2016). Understanding such interactions can be achieved through field observations, which may be severely limited in many cases. Simulation modeling has arisen in the last few decades as an efficient method to complement field investigations for understanding hydrologic processes. Simulation models are physically based hydrological models that simulate the terrestrial portion of the hydrologic cycle. The simulation output can determine discharge at the basin's outlets and along stream networks, hydrologic responses to land use and land cover changes, land–atmosphere interactions, water quality changes, aquifer response to groundwater extraction, and many other phenomena (Markstrom et al. 2008). These models calculate surface water and saturated/unsaturated subsurface water flow and other biogeochemical processes. A physically based hydrological model consists of mathematical descriptions of the surface process, subsurface process, external boundary conditions, internal boundary conditions, and initial conditions of a hydrologic system (Furman 2008). Physically based models have been applied in numerous studies to address a wide range of questions (Kampf & Burges 2007; Smith & Gupta 2012), such as climate change impacts/feedbacks (Loáiciga et al. 2000; Ferguson & Maxwell 2010; Sulis et al. 2011, 2012), stream-aquifer exchange (Peyrard et al. 2008; Frei et al. 2009), groundwater–lake interaction (Smerdon et al. 2007; Hunt et al. 2008), and agricultural sustainability (Schoups et al. 2005; Niswonger et al. 2017). Ebel et al. (2009) and Maxwell et al. (2014) provided a survey of studies that have applied hydrologic models for various purposes. There are numerous applications of hydrologic models, which may be classified based on peculiar characteristics and capabilities (Kollet et al. 2016). There are multiple categorizations of hydrologic models according to characteristics such as the number of hydrologic processes that are simulated, the types of mathematical conceptualization of hydrologic process, and the type of spatial representations of hydrologic elements and components (distributed, lumped, conceptual, etc.) (Barthel & Banzhaf 2016). This paper considers characteristics of hydrologic models related to their governing equations, the coupling of the governing equations, solution techniques, and numerical approaches. The governing equations for simulating the hydrological cycle are mostly based on mass and momentum conservation principles for surface and subsurface flow. The partial differential equations that are solved numerically express their complex nonlinear nature (Paniconi & Putti 2015). The governing equations of the hydrological models reviewed in this paper are based on the three-dimensional (3D) form of Richards' equation or combinations of the one-dimensional (1D) form of Richards' equation and a two-dimensional (2D) linear groundwater equation of variable saturated porous media that are coupled to some approximation of the Saint–Venant equations describing surface shallow flow in streams (Yen & Tsai 2001; Furman 2008).

Moreover, hydrological models can be classified into three different categories based on the surface coupling (SC) and surface governing (SG) equations as follows: (1) fully coupled, (2) iteratively coupled, and (3) loosely coupled (Panday & Huyakorn 2004; Furman 2008; Park et al. 2009; Rossman & Zlotnik 2013; Sebben et al. 2013; Maxwell et al. 2014). Fully coupled hydrologic models employ a numerical solution approach whereby, in every time step, the equations governing surface and subsurface processes are solved jointly in a transient simulation (Panday & Huyakorn 2004). Iteratively coupled hydrologic models feature governing equations of surface and subsurface processes that are solved separately. However, the cited equations establish feedbacks between them iteratively until numerical convergence is achieved. Loosely coupled hydrologic models feature surface and subsurface governing equations that are solved asynchronously such that outputs from one set of equation become inputs to another set of equations. Thus, each surface and subsurface equation is solved individually without receiving feedback from other equations (Guzha 2008).

There are three solution techniques for integrating surface and subsurface hydrologic processes: (1) first-order exchange (Panday & Huyakorn 2004), (2) continuity of pressure (Kollet & Maxwell 2006), and (3) the boundary condition switching procedure (Camporese et al. 2010). Equations governing hydrologic processes are generally based on systems of partial differential equations. Hydrologic models employ numerical solutions to solve those systems of equations. These numerical approaches can be the finite element method, the finite difference method, the finite volume method, or combinations of these methods.

There are a large number of hydrologic models, such as integrated hydrology model (InHM; VanderKwaak 1999), Catchment Hydrology (CATHY; Camporese et al. 2010), MIKE SHE (Graham & Refsgaard 2001), integrated groundwater-surface water model (IGSM; Labolle et al. 2003), HydroGeoSphere (HGS; Therrien et al. 2010), Parallel Flow (PARFLOW)-surface flow (Kollet & Maxwell 2006), WaterSHed Systems of 1D Stream-River Network, 2D Overland Regime, and 3D Subsurface Media (WASH123D; Yeh et al. 2004, 2006), GEOtop (Rigon et al. 2006), Penn State Integrated Hydrologic Model (PIHM; Qu & Duffy 2007; Kumar et al. 2009a, 2009b), Coupled Groundwater and Surface Water Flow Model (GSFLOW; Markstrom et al. 2008), FIHM (Kumar et al. 2009a, 2009b), IRENE (Spanoudaki et al. 2009), Cast3M (Weill et al. 2009), PAWS (Shen & Phanikumar 2010), and Open-GeoSys (Delfs et al. 2009; Kolditz et al. 2012) to name a few. A few commonly used models are briefly reviewed in the following subsections. Tables 1 and 2 summarize the main features of the above-mentioned hydrologic models, which are commonly used.

The MIKE SHE model was developed under the name Système Hydrologique Européen (SHE) in 1982. It is a proprietary model. This model was sponsored and developed by three European organizations: the British Institute of Hydrology, the French consulting company SOGREAH, and the Danish Hydraulic Institute (DHI) (Taghavi et al. 2013). MIKE SHE is a distributed physically based model introduced in 1990. MIKE SHE is an integrated hydrological modeling system with a structure covering a wide range of hydrologic processes such as evapotranspiration, overland flow, unsaturated flow, groundwater flow, and channel flow, designed for water managers. The hydrological processes are represented at various levels of complexity in time and space. MIKE SHE may be linked to the river hydraulics model MIKE 11 and ECOLAB model for water quality simulation. MIKE SHE features a Windows-based user interface that is used for water quality, water budget analysis, and parameter estimation. It provides users with calculation of water demands based on crop and determines water allocation through the irrigation decision support system. MIKE SHE describes fluid motion by transient partial differential equations that are solved numerically (Demetriou & Punthakey 1999; Shi et al. 2012). The MIKE SHE's code is found in the user's guide.

PIHM is described as a fully coupled hydrologic model with two multi-scale and multiprocessing features. PIHM uses a semi-discrete finite volume approach for coupling groundwater, land surface components, and surface water to make fully coupled ordinary differential equations from flow systems described by partial differential equations (Qu & Duffy 2007). PIHM simulates processes such as infiltration, evaporation, transpiration, streamflow, subsurface flow, overland flow, and recharge (Seo et al. 2018). The PIHM system is open-source software.

This finite element-based model was developed in the 1970s and was expanded in the 1990s with its application to California's Central Valley GW-SW Model in support of the Central Valley Project Improvement Act (Labolle et al. 2003). IGSM is an integrated model that simulates the complete hydrologic cycle within a basin. The IGSM has been used in various projects undergoing code modifications for new features (WRIME 2003). There is a high similarity between the IGSM and Streamflow Routing Package of MODFLOW. The model solves confined and unconfined flow equations.

GEOtop is a distributed hydrological model that simulates 3D coupled water and energy (heat) budgets on and beneath the soil surface, considering both turbulent and radiative fluxes within a grid (Endrizzi et al. 2014). GEOtop applies a numerical solution of the 3D Richards' equation to model saturated and unsaturated subsurface flow. The kinetic wave approximation offered by Gottardi & Venutelli (1993) describes the surface water flow. Rigon et al. (2006) showed that GEOtop might be easily interfaced with climate and bounded-area meteorological models to simulate river basin hydrology.

The HGS model was developed by Therrien et al. (2010). HGS is a 3D model with fully coupled surface–subsurface simulators. It fully couples surface flow, water quality, and subsurface flow processes. HGS applies a 2D diffusive wave approximation of the Saint–Venant equations for surface flow, the 3D Richards' equation for unsaturated/saturated flow simulation, and a classic advection-dispersion equation for solute and thermal energy transport. Besides, hydraulic features such as streams, rivers, subsurface wells, water supply lines, and drain pipes are simulated by 1D analytical and empirical equations, such as the Hagen–Poiseuille, Manning's, and the Hazen–Williams formulas. HGS's numerical solver implements a control-volume finite element method. HGS can apply multiple spatial discretization options ranging from simple rectangular domains to irregular domains featuring intricate geometry and layering. HGS considers depression storage such as rill storage and storage exclusion for rural and urban environments. There are no artificial storage control features in HGS at present (Aquanty 2013). HGS applies adaptive time-stepping, whereby the simulation time step may change through algorithmic execution guided by the rate of change of the dependent variables. HGS applies a Newton–Raphson linearization method to handle nonlinear governing equations. HGS features modules that calculate evaporation from bare soil and water bodies, vegetation-dependent transpiration with root uptake, snowmelt, and soil freeze/thaw. The computational efficiency of HGS relies on a parallel computational framework applying OpenMP (Hwang et al. 2014). HGS's code is written in FORTRAN 95 and is not open source. It can be used on any Microsoft Windows or Linux based platforms. HGS applies a graphical user interface (GUI) for grid generation and HSPLOT and Tecplot for 3D virtualization and animation. In addition, HGS has interfaces with GIS tools such as ArcView and ArcInfo for enhanced spatial data analysis.

PARFLOW is a fully integrated and parallel watershed flow model introduced by Kollet & Maxwell (2006). PARFLOW integrates a land surface model (the Common Land Model (CLM)) and a groundwater model (PARFLOW). The model features a fully coupled 2D kinematic wave approximation of the Saint–Venant equations for overland flow and the 3D form of Richards' equation. The pressure continuity approach is implemented for the solution of the coupled surface and subsurface equations. CLM processes consider biogeochemistry, land energy budget, and snow processes. CLM computes the parameters as follows: momentum, sensible heat, latent heat, ground heat fluxes, outgoing long-wave radiation, and the surface albedo. PARFLOW handles the nonlinearity of governing equations with the Newton–Krylov approach of multigrid preconditioning (Maxwell et al. 2014). PARFLOW can be run on a range of serial or parallel Linux, Unix, and OSX platforms. PARFLOW was originally written in C and updated to FORTRAN 90/95. PARFLOW has no visualization capabilities at present. The codes are open source. PARFLOW does not contain any lake or reservoir storage modules (Maxwell et al. 2009).

GSFLOW was developed by the United States Geological Survey (USGS; Markstrom et al. 2008). GSFLOW couples the 2D precipitation-runoff modeling system (PRMS) with the 3D modular groundwater flow model MODFLOW-NWT. The latter model is an upgrade of the 3D groundwater flow model MODFLOW-2005, which consists of several packages to simulate unsaturated zone flow (UZF1) (Niswonger et al. 2017; Niswonger 2020), Streamflow Routing (SFR1) (Prudic et al. 2004), lake–aquifer interactions (LAK3) (Merritt & Konikow 2000), groundwater extraction and artificial recharge, unsaturated flow precipitation recharge, subsurface drains, and evapotranspiration. GSFLOW assumes vertical 1D unsaturated flow in the vadose zone coupled to 3D groundwater flow. The PRMS and MODFLOW-NWT are coupled by an iterative Newtonian linearization method. PRMS simulates parameters such as runoff, evapotranspiration, interflow, and infiltration by balancing mass and energy budgets of the snowpack, plant canopy, and soil zone based on distributed climatic variables, which the most important ones are as follows: (1) temperature (2) precipitation, and (3) solar radiation. PRMS applies an approximation of kinematic wave with single-direction and unsteady surface flow (Leavesley et al. 2006). In addition, PRMS approximates the effects of the urbanization on groundwater recharge. GSFLOW calculates all major components of the hydrologic cycle in the saturated and unsaturated zones, streams, and lakes. GSFLOW does not contain modules for handling hydraulic storage structures. It features a dynamic control option for the diversion of streams. Outflow from the lake package can be specified overtime or controlled with a weir or spillway.

GSFLOW was written in Fortran 90 and C languages, and the codes can be run in Microsoft Windows and Linux platform. GSFLOW codes are open source and can be run in three modes, namely, the GSFLOW mode, only PRMS mode, and only MODFLOW mode. GSFLOW runs on a daily time step (Markstrom et al. 2008).

CATHY is a physically based catchment-scale model developed by Camporese et al. (2010). The model is based on coupling the subsurface module FLOW 3D and the surface module SURF-ROUT. The former module is a 3D form of Richards' equation solved by finite elements for variably saturated subsurface flow. The latter module is a digital elevation model (DEM)-based surface runoff module to simulate slope and stream transport, including lakes and pools. Surface flow is routed on a path-based 1D diffusion wave equation and the Muskingum–Cunge routing technique solved by finite differences. The coupled solution of the surface and subsurface equations is accomplished with a boundary-condition switching process, which simulates physical processes by switching boundary conditions automatically between atmospheric-controlled fluxes to soil-controlled (constant head) conditions through infiltration excess and saturation excess while calculating changes in surface water storage.

MODHMS was developed by Panday & Huyakorn (2004). It is a spatially distributed model that simulates surface and subsurface flow, and flows through a network of channels that takes all relevant processes and hydrologic interactions into account. The model is based on coupling the 3D form of Richards' equation for variably saturated flow, the 2D diffusion wave approximation for overland flow, and the 1D diffusion wave approximation for surface water features such as rivers, channels, pipes, lakes, reservoirs, and ponds. Surface and subsurface evapotranspiration is calculated based on physically based formulations. The model simulates hydraulic structures such as weirs, culverts, manholes, drop-structures, bridges, gates, and pumps within the channel flow regime. A set of operational rules can manage hydraulic structure. MODHMS relies on the first-order exchange method for the fully coupled solution of the governing equations. The subsurface equations are solved with the finite difference method. The channel flow equation is solved with the finite volume method.

The InHM was introduced by VanderKwaak (1999). It is a physically based model that fully couples surface flow and variably saturated porous media with or without macrospores and fractures to simulate surface and subsurface flows. InHM considers streamflow generation in addition to dissolved chemical transport (solute transport), but it does not include snowmelt, canopy interception, and sediment transport processes. InHM applies a 2D diffusion wave approximation of the Saint–Venant equation to simulate overland flow and channel flow. It implements a 3D form of Richards' equation fully coupled to surface flow. Linkage of surface and subsurface flow is based on the first-order exchange coefficient method.

WASH123D is an integrated numerical watershed model developed by Yeh et al. (2004, 2006). WASH123D can simulate multiple processes as follows: (1) hydrologic fluxes such as evapotranspiration, evaporation, infiltration, recharge; (2) fluid flow such as land surface runoff, river/stream/canal networks hydraulics and hydrodynamics, vadose zones' interflow and saturated zones' groundwater; (3) salinity and thermal transport (in both surface water and groundwater); (4) surface water's sediment transport; (5) water pollutants transport (any number of reactive constituents); (6) biogeochemical cycles include nitrogen, phosphorous, carbon, and oxygen; and (7) biota kinetics such as phytoplankton, zooplankton, algae, coliform, plants, and bacteria. Flow through streams, canals, open channel networks, and overland flow are simulated with the 1D and 2D of Saint–Venant equations employing total kinematic wave, diffusive wave, and dynamic wave approximation options, respectively. Flow through saturated and unsaturated porous media is calculated with the 3D form of Richards' equation. Salinity, thermal, and sediment transport in river networks and overland regime are simulated with modified advection-dispersion equations featuring empirical formulas for erosion and deposition. Solute transport is simulated by means of the advection-dispersion-reaction equations with reaction-based mechanistic approaches to water quality modeling. The numerical method applied to approximate the kinematic wave, diffusive wave, and dynamic wave, and Richards' equations are semi-Lagrangian, Galerkin finite element or the semi-Lagrangian, and Lagrangian–Eulerian finite element, and Galerkin finite element methods, respectively. The length scale of modeling ranges from meters to thousands of kilometers for small dam-break problems to large watershed simulations, respectively. The time scale of modeling ranges from seconds (for dam-break problems) to tens of years (for simulation of large watersheds). WASH123D includes flow control structures such as weirs, gates, culverts, pumps, levees, and storage ponds featuring operational rules for pumps and control structures. The WASH123D code consists of eight modules to cope with multiple media. WASH123D is coded in Fortran77, Fortran90, and runs in Linux and Windows platforms. Reviews of other widely applied hydrologic models (e.g., Hydrologic Simulation Program Fortran (HSPF), the Soil Water and Assessment Tool (SWAT), Hydrologic Engineering Center Hydrologic Modeling System (HEC HMS), the stormwater management model (SWMM), integrated modeling platforms featuring data management and input modules, hydrologic, water quality, and sediment transport models, and data analysis and visualization modules (the Better Assessment Science Integrating Point and Non-point Sources (BASINS), among others) can be found in ASCE (2017).

Koch et al. (2016) compared the HGS, MIKE SHE, and PARFLOW-CLM models' capabilities to simulate spatiotemporal discharge and soil moisture in a small catchment in Germany. They showed that the 3D capability for simulating flow in the unsaturated zone did not essentially improve the predictions of soil moisture in the study catchment. Maxwell et al. (2014) compared the ability of seven hydrologic models (CATHY, HGS, OGS, PIHM, PARFLOW, PAWS, and RIBS-VEGGIE) using benchmark problems. Their results established good inter-model agreement of different models' results for simpler benchmarks (runoff generation dynamics). Nevertheless, the agreement of results decreased as benchmarks became more complex, beset by heterogeneity and phreatic-surface dynamics. Sulis et al. (2010) compared two hydrologic models (PARFLOW and CATHY) to simulate GW-SW interactions. Results showed a good agreement on hydrograph shape and flow distribution with depth calculated with the two models.

Hydrologic models have various capabilities to simulate physical behavior of surface water and groundwater (Kollet & Maxwell 2008; Ferguson & Maxwell 2010; Rozemeijer et al. 2010; Shen et al. 2013; Singh 2014; Rajagopal et al. 2015; Mateus & Tullos 2016). However, these models cannot determine best management scenarios. In other words, simulation models are useful in answering ‘what if’ questions and simulating answers for management problems like: what is likely to happen if a particular operating policy is implemented? Hydrologic models evaluate the physical behavior of a hydrologic system under specified conditions. Yet, hydrologic simulations are not well suited or designed for identifying best management solutions for a specific goal (Loucks et al. 1981). Rather, coupling optimization methods with hydrological flow models seem to be necessary to achieve optimal management policies. The optimization algorithms adjust the model's parameters to improve an objective function, thus yielding optimal solutions to conjunctive management problems.

Optimization methods have been widely used in water resources systems to find optimal management strategies for various purposes such as controlling non-point agricultural contamination, determining optimal location of low impact development (LID) infrastructures for urban runoff control, enhancing the yield via water demand optimization, and preventing seawater intrusion by optimizing withdrawal in coastal aquifers (Cunha & Sousa 1999; Koech et al. 2014; Aboutalebi et al. 2016; Amirkhani et al. 2016; Deng et al. 2016; Garousi-Nejad et al. 2016; Soleimani et al. 2016a, 2016b). Optimization methods are classified as gradient-based methods and evolutionary or metaheuristics methods. The former methods include linear programming (LP), nonlinear programming (NLP), and dynamic programming (DP). The latter methods include hundreds of methods, among which the pioneering GA (Holland 1975) is widely known, and do not require derivative calculations. However, the large number of iterations makes the optimization process time-consuming, particularly in combination with complex GW-SW models. That is why surrogate modeling has been introduced. A summary of optimization methods and their applications in the conjunctive use of GW-SW resources follows.

LP is a relatively simple technique to optimize a linear objective function subject to linear equality and inequality constraints. Most real-world water resources problems are nonlinear. Applying LP to real water recourse problems requires the introduction of simplifications that reduces the accuracy and applicability of results (Bertsimas & Tsitsiklis 1997). Nevertheless, there are multiple applications of LP to model stream-aquifer interactions and the conjunctive use of GW-SW resources (Morel-Seytoux 1975; Tyagi & Narayana 1981; Tyagi & Narayana 1984; Gaur et al. 2011; Lu et al. 2011).

LP has limitations in solving GW-SW problems (Gorelick et al. 1984). In contrast, NLP can tackle many nonlinear problems that are amenable to mathematical differentiation and gradient calculations in convex or concave objective function spaces. Time-consuming and complex calculation of derivatives (gradients) may become too cumbersome; in some cases, this is a key limitation of the NLP method (Wismer & Chattergy 1978). Another limitation occurs when the GW-SW problem features non-differentiable functions or discrete variables. There have been many applications of NLP to the conjunctive use of GW-SW (Haimes & Dreizin 1977; Khan 1982; Ghahraman & Sepaskhah 2004; Chiu et al. 2009; Montazar et al. 2010; Huang et al. 2012).

One classical optimization method is DP, developed by Bellman (1957). DP breaks down a multiple-decision-making problem into a sequence of simpler subproblems relying on the principle of optimality. By solving the subproblems, one arrives at the solution of the entire optimization problem. The DP sequential search algorithm can solve nonlinear, non-differentiable, discontinuous, discrete variable, deterministic or stochastic problems. As the number of states and decision variables increases, the computational burden grows exponentially, leading to the so-called curse of dimensionality. GW-SW problems involving continuous variables must be discretized by approximating the continuous variables with a set of discrete variables, which compounds the computational burden and reduces the accuracy of the solutions. Nevertheless, the advent of powerful computational technology has allowed many applications of DP to the solution of GW-SW resources problems (Hall et al. 1968; Burt 1970; Rao et al. 1992; Umamahesh & Sreenivasulu 1997; Philbrick & Kitanidis 1998; Naadimuthu et al. 1999; Ben Alaya et al. 2003; Karamouz et al. 2004; Tran et al. 2011). Numerous studies applied dynamic programming in conjunctive use planning on account of its ability in undertaking sequential decision-making and incorporating stochasticity of hydrological procedures.

Evolutionary and metaheuristic algorithms are based on computational and search methods influenced by artificial intelligence (AI). A heuristic function, also called simply a heuristic, is a function that ranks alternatives with search algorithms at each branching step based on available information to decide which branch to follow (Newell & Simon 1976; Bozorg-Haddad et al. 2017). The GA, particle swarm optimization, and shuffled complex evolution are examples of evolutionary and metaheuristic algorithms. The GA is the pioneering evolutionary optimization algorithm introduced by Holland in 1975. They reach near globally optimal solutions with acceptable errors in problems that some classical methods cannot solve, such as nonlinear, non-differentiable, discontinuous, mixed discrete/real variables, non-convex, deterministic, stochastic, single objective, and multi-objective formulations. Evolutionary and metaheuristic algorithms are less likely to converge to local optima than classical methods (gradient-based NLP methods). However, they require numerous evaluations of the objective function, which may introduce a heavy computational burden solving coupled physically based simulation models and evolutionary algorithms. In addition, most evolutionary algorithms require the specification of parameters that may require calibration.

The optimization with evolutionary and metaheuristics algorithms begins with a randomly generated population of management policies. The water resources system is simulated with hydrologic–economic models based on the current population of generated management policies. The performance of the water resources system is evaluated with single objective or multi-objective functions. The evolutionary algorithm then creates an improved population of management policies with which to simulate anew the performance of the water resources system. This iterative process continues until the performance of the water resources system cannot be improved any further. At that point, the simulation–optimization method has identified an optimal management policy and the corresponding response of the water resources system.

This modeling approach converts complex functions with much simpler ones in an iterative model assessment procedure. The surrogate modeling approach is divided into two types. The first type discovers the relationship between multiple descriptive variables and a model output variable. Support Vector Machine, Kriging, Probabilistic Collocation Method, Radial Basis Function, and Dynamic Coordinate search using Response Surface models are typical samples of the first type. The second type of surrogate modeling reduces the order of the original complex model. The second type has been coupled to surface water models and groundwater models separately; yet, its application in fully integrated SW-GW models is rare (Wu et al. 2016). A few investigations have combined surrogate modeling methods for the purpose of enhancing the accuracy of results. Numerous studies have used single or multiple surrogate models and have compared the computational costs with those of direct optimization, which in some cases has produced up to 90% reduction in the computational burden (Christelis et al. 2019).

Simulation–optimization methods provide an in-depth understanding of the effect of management strategies on water resources systems. They systematically evaluate management/operational options, including the optimal choices via search algorithms. Nevertheless, coupling complex GW-SW models with optimization tools is challenging because the derivatives of the functions that appear in the coupled models are not determined analytically. Also, optimization with evolutionary algorithms is time-consuming. Therefore, except for a few studies incorporating GW-SW models in optimization, most published works apply surrogate modeling. Simulation models and optimization tools are coupled through embedding techniques or response matrices. The embedding method discretizes the flow equations as constraints into a linear program, and the simulation model is solved within the optimization process. The response matrix method, on the other hand, performs simulation and separation in two separate steps. Firstly, the response of the aquatic system to hydraulic stress such as groundwater withdrawal is computed by the external simulation model. Next, the response matrix is incorporated by applying superposition and linear system theory (Singh 2014). Alternatively, nonlinear optimization can be applied by means of evolutionary and metaheuristic algorithms, dynamic programing, or NLP as deemed appropriate. Examples of studies that have coupled physically based models to optimization models to manage surface and groundwater resources are listed in Table 3.

In addition to the works listed in Table 3, there are numerous studies indicating the benefits of optimization in conjunctive water use planning, while applying other simulation models than those mentioned in the ‘Hydrologic simulation models’ section. Brookfield & Gnau (2016) coupled the HGS model and a surface water operation model called Operational Analysis and Simulation of Integrated Systems (OASIS) to evaluate various water management strategies, such as structural and operational strategies, for meeting future demands of a water-stressed agricultural basin in the USA's Lower Republican River. OASIS applied the LP method to minimize agricultural water demand deficit.

Insofar as the application of evolutionary algorithms is concerned, Karamouz et al. (2004) proposed a simulation-based DP optimization model for conjunctive GW-SW planning and management in Iran that minimizes irrigation water supply shortages and pumping costs and controls the average groundwater table fluctuations. The optimal time schedule of water diversion was calculated for areas dependent on groundwater. Yang et al. (2009) presented a model that coupled a groundwater simulation model and DP for conjunctive, multi-objective, GW-SW management in Taiwan that minimized fixed and operating costs. Safavi et al. (2010) applied a trained artificial neural network (ANN) model to simulate the interaction between surface water and groundwater linked to the GA as the optimization model. The ANN-GA coupled algorithm minimizes shortages of irrigation water subjected to constraints on the cumulative groundwater drawdown and on the maximum capacity of surface irrigation systems. Chang et al. (2013) developed a fuzzy inference system (FIS) for managing conjunctive GW-SW. Safavi & Rezaei (2015) combined the FIS and fuzzy neural networks (FNNs) for simulation with the multi-objective GA for optimization in a simulation–optimization for surface water and groundwater resources in the Najafabad plain in central Iran. Heydari et al. (2016) coupled the ANN model for groundwater level simulation and the genetic programming model for TDS concentration prediction with the multi-objective non-dominated sorting GA NSGA-II for modeling water allocation with the simulation–optimization modeling approach. Rezaei et al. (2017) employed fuzzy multi-objective particle swarm optimization (f-MOPSO) to the management of GW-SW. f-MOPSO is relatively simple in concept, easy to implement, and computationally efficient for solving large-scale optimization problems.

Most published studies have solved the conjunctive use management problems, applied classical optimization methods and reduced the computational burden by simplifying the groundwater equations (Karamouz et al. 2004). Safavi & Esmikhani (2013) proposed a simulation–optimization model for analyzing the conjunctive use of surface water and groundwater in the Zayanderud basin, Iran. They applied a groundwater model (MODFLOW) and obtained groundwater levels; this was followed by applying a data mining method (support vector machine, SVM) as a surrogate of the MODFLOW model. The SVM model was linked to the GA to extract optimal surface water releases and groundwater withdrawal to meet agricultural water demand. Wu et al. (2015, 2016) reported methodology to optimize complex conjunctive use problems. The latter authors proposed a surrogate-based optimization for integrated surface water-groundwater modeling, which constitutes an integrated model for optimizing conjunctive river-aquifer management.

The water resources system literature indicates the importance of searching for reservoir operational rules and groundwater extraction plans to optimize conjunctive GW-SW resources management. Optimal management of GW-SW resources, particularly in areas dependent on groundwater resources, must rely on coupled simulation–optimization linking hydrologic models with optimization algorithms that accurately represent real-world conditions. A review of coupled physical-based hydrologic models and optimization models to achieve optimal conjunctive management GW-SW resources is presented in this paper. This paper's review of the most used GW-SW models, their features, and suitable optimization methods provides valuable clues for selecting among them for the purpose of research, planning, and policymaking. The application of simulation–optimization has clear benefits for planners, such as determining: (1) the optimal conjunctive management of GW-SW resources, (2) the best possible management of reservoirs and groundwater withdrawal during droughts, (3) the optimal withdrawal to prevent seawater intrusion in coastal aquifers, (4) the indices of reliability, resiliency, and vulnerability that optimize the conjunctive management of GW-SW resources, and (5) the reservoir releases that optimize hydropower production as well as groundwater and streamflow interactions downstream of reservoirs. There are few published articles in the aforementioned categories, and new optimization methods are not commonly applied. Other fields in water resource management could benefit from simulation–optimization tools. The following fields of inquiry are herein identified as priorities for future research:

• Most previous studies have focused on surrogate-based modeling and evolutionary algorithms; however, application of gradient-based methods, coupling the GW-SW simulation models with different optimizers, and comparing the results of various optimizers has not been adequately covered. Investigation on this subject would provide helpful insights for watershed planers.

• The existing literature is predominantly abundant with coupled GW-SW simulation–optimization tools focused on irrigation management, hydropower production planning, and coastal aquifer management, with less attention given to stormwater management considering groundwater quality impacts. GW-SW simulation–optimization tools can be applied to optimize stormwater control strategies and assess water quality and quantity effects on the conjunctive use of surface water and groundwater.

The authors thank the Iran National Science Foundation (INSF) for its support of this research.

None.

The data that support the findings of this study are available from the corresponding author upon reasonable request.