\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Optimization of Groundwater Pumping and River-Aquifer Exchanges for
Management of Water Resources}



\author[1]{Mayank Bajpai}%
\author[1]{Shishir Gaur}%
\author[1]{Anurag Ohri}%
\author[1]{Shreyansh Mishra}%
\author[2]{Hervé Piégay}%
\author[3]{Didier Graillot}%
\affil[1]{Indian Institute of Technology BHU Varanasi}%
\affil[2]{CNRS}%
\affil[3]{Ecole des Mines de Saint-Etienne}%


\vspace{-1em}



  
  \date{June 02, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Groundwater pumping influences the rate of River-Aquifer (R-A) exchanges
and alters the water budget of the aquifer. Therefore, fulfilling the
total water demand of the area, with an optimal pumping rate of wells
and optimal R-A exchanges rate, is important for the sustainable
management of water resources and aquatic ecosystems. Meanwhile,
comparison of the output of different simulation-optimization
techniques, which is used for the solution of water resource management
problems, is a very challenging task where different Pareto fronts are
compared to identify the best results. In the present work, mathematical
models were developed to simulate the R-A exchanges for the lower part
of the River Ain, France. The developed models were coupled with
optimization models in MATLAB environment and were executed to solve the
multi-objective optimization problem based on the maximization of
pumping rates of wells and maximization of groundwater input into the
river Ain through R-A exchanges. The Pareto front developed by different
simulation-optimization models was compared and analyzed. The Pareto
fronts were juxtaposed based on the convergence, total diversity, and
uniformity with the help of different performance metrics like
hypervolume, generational distance, inverted generational distance, etc.
The impact of different groundwater models based on domain size and
boundary conditions was also examined. Results show the dominance of
MOPSO over other optimization algorithms and concluded that the
maximization of pumping rates significantly changes after considering
the R-A exchanges-based objective function. It is observed that the
model domain also alters the output of simulation-optimization,
therefore the model domain and corresponding boundary conditions should
be selected carefully for the field application of management models.
ANN models were also developed to deal with the computationally
expensive simulation model by reducing the processing time and were
found efficient. Keywords: Simulation-Optimization, Multi-Objective
optimization, Artificial Neural Network, River-Aquifer exchanges.%
\end{abstract}%



\sloppy


\textbf{1. Introduction:}

River-aquifer (R-A) exchange influences both quantity and quality
aspects of groundwater and river water systems significantly.
Consequently, proper and effective quantification and representation of
river-aquifer exchange are very important for the management of water
resources and aquatic ecosystems. Climate and groundwater pumping for
irrigation has caused rapid groundwater depletion in India (Dangar et
al.,2021). Several works have been carried out about the identification
and the quantification of groundwater and surface water interactions
i.e. R-A exchanges (Constanzt, 1998; Sophocleous, 2002; Becker et al.,
2004; Anderson, 2005; Kalbus et al., 2006; Keery et al., 2007; Lowry et
al., 2007; Yang et al., 2017). Some studies on R-A exchanges also
demonstrate the effect of it on river temperature (Westhoff et al.,
2007; Burkholder et al., 2008; Hebert et al., 2011) along with the
emphasis on exchanges in hyporheic zone (Burkholder et al., 2008).

Simulation-optimization (S-O) techniques are often used in identifying
the optimal management practices of groundwater for the selected region.
Many challenges are associated with this S-O technique like the
selection for proper simulation or optimization algorithm and comparison
of the results of various techniques used. In the multi-objective
optimization problems, comparison of results becomes more significant as
Pareto fronts generated in this approach can have some similarities with
some diversity. Jha et al. (2020) statistically evaluated the
relationship between groundwater pumping rates and groundwater levels
during pre-monsoon and post-monsoon seasons using a numerical
groundwater-flow simulation model. The simulation-optimization (S-O)
model used linear programming (LP) optimization to handle groundwater
hydraulic optimization management problems after convolution techniques
were employed to integrate the simulation model with an optimization
method. Kamali Asghar (2017) used a mathematical simulation optimization
programming model to handle the issue of aquifer management, relying on
the stability of water quality and quantity. The simulation-optimization
models and the construction of the sustainability index for determining
the best point of the Pareto front were used in the study's modeling.
Abd-Elmaboud et al. (2021) introduced a new model that connects
geomorphological and hydrogeological data with recharge rates. The
MFUSG-PSO model with an indirect simulation-optimization technique for
the inversion of the groundwater flow problem was used to calibrate the
recharge rates. The second phase involved training a CFNN model to
connect the calibrated recharge rates with freely available
geomorphological and hydrogeological data to construct
interrelationships. For the best possible conservation of water
resources, Conant et al. (2019) suggested that it is critical to
understand and quantify exchange activities between rivers and
groundwater. The exchange~between rivers and groundwater is significant
in a variety of current concerns, including providing drinking water,
characterizing and managing environmental flow regimes, preserving or
restoring riverine ecosystem health and functioning, and alleviating
toxins. For reaches in plains with flow monitoring data. Li et al.
(2020) introduced a cumulative exchange fluxes method based on surface
water balance to study GW-SW interactions. The dynamic change processes
of GW-SW interactions can be qualitatively and quantitatively judged
through a curve of cumulative exchange fluxes by this method. The
hyporheic zone, which is usually regarded as a biogeochemical hotspot,
is where surface-groundwater interactions are most frequent and
according to McClain et al. (2003), these hotspots are the key to manage
water resources effectively.

Calandra et al. (2014) proposed an approach to make full use of the
simulation algorithms developed by Bayesian Optimization to analyze the
quality of Pareto fronts, unlike the existing algorithms which only
returned the Pareto solution sets without considering the qualitative
assessment. The approach was able to delineate the actual Pareto fronts
better in presence of noise and at the same time can also perform
sensitivity analysis of the parameters concerning the quality of model
output. Lobato et al. (2016) have presented a Bayesian method called
PESMO to identify Pareto solution set to Multi-Objective optimization
problems. The evaluation points were chosen as such to minimize the
entropy of the Pareto set. When compared with existing techniques, PESMO
was found to produce better results with a smaller count of evaluations,
while the decoupled estimation led to an improved performance
particularly when the existing techniques lose efficiency with an
increase in objectives. Emmerich et al. (2018) have also discussed the
fundamentals and evolutions in the field of Multi-objective
optimization. The topics covered in their article include-
order-theoretical foundations, scalarization approaches, and optimality
conditions. In context to the evolutionary methods, three
state-of-the-art techniques were discussed namely NSGA-II, SMS-EMOA, and
MOEA/D. NSGA-II representing the Pareto-based approach, SMS-EMOA
exemplifies the Indicator-based approach while MOEA/D is an example of a
decomposition-based approach. The choice of the correct technique
depends on the number of objectives, count of solution sets, desired
distribution of the solutions along with the location and shape of the
Pareto front. Belakaria et al. (2019) have addressed the optimization of
the Multi-Object (MO) black box to determine true solutions of
Pareto-set by reducing the count of function evaluations. The study
proposed a new algorithm called Max-value Entropy Search for
Multi-objective Optimization (MESMO) to ascertain an optimized design
that efficiently trades off among performance, power requirement, and
area overhead. Their approach used the output matrix entropy function to
efficiently select the input parameters, to obtain highly accurate
Pareto-set solutions. The algorithm was found to be constantly
outperforming the existing state-of-the-art algorithms in computing
Pareto set solutions.

Binois et al. (2014) have taken into account Kriging metamodel for
estimating the Pareto front and also to quantify bias associated with
the solution set at any phase of the multi-objective optimization. The
approach taken by them assumed the original dataset to be having
Gaussian distribution. The concept of random set theory has been used to
compute Vorob'ev deviation to capture the variability of the dataset
with non-dominated points. This method applied on several numerical
problems yielded satisfactory output in accurately determining of Pareto
front. Horn et al. (2016) proposed a model based on Random Forests for
model-based multi-objective optimization to determine the Pareto front
for mixed hierarchical configuration problems. A two-phase parameter
experiment was carried out. The results of the single model experiment
were quite promising, however, the solution derived through a
multi-model approach could be improved further. Especially, the bias
estimation associated with the Random Forest model needed to be
resolved.

Cao et al. (2017) have tried to quantify the bias of the estimated
solution of multi-objective problems using two versions of the
normalized hypervolume. A case study was carried out in a chemical
process with three sets of solutions, each having a different
optimization goal. Results showed that the normalized hypervolume showed
great accordance to the ideal situations while taking into account the
degree of the trade-offs between the possibilities considered to compute
the front. Bassi et al. (2018) have considered a new approach to
evaluate statistical parameters for a set of objects equivalent to
surfaces and arcs. Their approach is based on the identification of the
most representative member of a family tree and was found useful to
address the uncertainties associated with Pareto fronts. The Pareto
front was ascertained by minimizing the hypervolume between the front
and the reference points. The algorithm was tested on a complex 5 bar
truss structure with satisfactory outputs being obtained. Avent et al.
(2019) proposed a new approach called DPARETO, to determine the
trade-off among various differentially cloistered algorithms. Bayesian
Optimization (BO) was used to concurrently optimize both privacy and
utility parameters of a Pareto front. Moreover, they also established
the effectiveness of BO in creating visualization interfaces helpful in
decision making. Bionis et al. (2019) have designed a GPareto library
for R which enables the optimization of multi-objective algorithms for
functions associated with a black box. Moreover, the package also
contained several algorithms for the accurate assessment of the
associated bias. The study also proposed several infill criteria in
minimizing the bias associated with several optimization models such as
efficient global optimization technique.

Marjit R. (2009) has made a noble attempt to increase the efficiency of
the Building performance simulation (BPS) tool by combining it with a
robust multi-objective optimization algorithm. The thesis aimed at
developing a multi-objective algorithm with a meta-model to optimize the
simulation of a Pareto front. Asadzadeh et al. (2014) on the other hand,
have developed a new technique called Convex Hull Contribution (CHC) to
address multi-objective (MO) optimization functions associated with
Pareto fronts. They demonstrated the effectiveness of CHC in enhancing
the accuracy of Pareto archived multi dimensioned search while
addressing the multi-objective problems. Audet et al. (2018) have
presented a review of the algorithms devised so far to address the
multi-objective-based optimization in the determination of Pareto
fronts. An assessment of a large number of performance indicators was
carried through various algorithms. Total 57 performance indicators were
grouped into four categories by their properties, viz. cardinality,
convergence, distribution, and spread. Hollermann et al. (2019) have
proposed a new flexible method to automatically identify one optimum
design for multi-objective decision-making. Their approach addresses
both economic and environmental aspects to identify a viable design. The
bias associated with the parameters while predicting a future event
could be handled easily by developing an extension of this algorithm.

Present work was carried out to address the challenging issue of R-A
exchanges which can be analyzed through numerical models. The output of
the numerical models depends on the model domain considered, different
hydrogeological parameters used and boundary conditions applied.
Considering the R-A exchanges in groundwater resource management
problems can help in finding out more efficient and river health
inclusive management practices through simulation-optimization models.
Meanwhile, selection of the optimization algorithm, model domain
demarcation and the comparison of optimum results to provide a
physically meaningful management strategy is still a bottleneck in the
groundwater management problems. Therefore, different model domain
demarcations for the same river system were considered for the
development of the groundwater model, and their impact on R-A exchanges
was analyzed. In addition, different optimization techniques were also
used and their outputs i.e. Pareto fronts were compared. Finally, this
set of optimal results were interpreted in terms of groundwater
management.

\textbf{2. Study Area:}

The Ain River has an oceanic hydrologic regime and drains a basin of
about 3630 km\textsuperscript{2}. It is located in the southern Jura
Mountains, France. The length of the Ain River is about 200 km and it is
the right bank tributary of the Rh\selectlanguage{ngerman}ône River. The lower part of the Ain
River, which is located between the Allement dam and the confluence with
the Rhône River, was taken for the study purpose. Figure 1 shows the
location of the lower Ain river area. Along this section, the average
channel width is about 60 m and the slope is about 1.3\%. The mean
annual discharge of the river, between 1959 to 2019, was evaluated as
120 m\textsuperscript{3}/s at the measuring station of Chazey-sur-Ain.
Most floods occur between October and March (probability of occurring
flood discharge of 1200 m\textsuperscript{3}/s, once in 10 years, at
Chazey-sur-Ain) and summer is characterized by low flow conditions.
Further detail about the study area is provided in the groundwater model
section.

\textbf{3. Methodology}

Two major objectives were addressed in this study. The first objective
was to identify the maximum pumping rates of the wells by considering
the R-A exchanges component. This also includes examining the impact of
different types of model domains on the output of R-A exchanges and to
find the best model domain and boundary conditions. Whereas the second
objective was to compare the Pareto front developed by different coupled
simulation-optimization models. Different groundwater model domains were
taken by using the different boundary conditions and the impact of each
domain on R-A exchanges and correspondingly the water budgeting of the
aquifer was examined. Figure 2 shows the flow chart of the methodology.
Total four optimization techniques were compared i.e. Multi-objective
Genetic Algorithm (MoGA), Multi-objective Particle Swarm Optimization
(MOPSO), Pareto Search (PS), and Multi-Objective Evolutionary Algorithm
based on Decomposition (MOEA/D). To deal with a large number of wells,
the wells in the domain were grouped based on the municipal zones and
distance of the wells from the river. If any well in the municipal zone
had a distance from the river less than a threshold of 1 km, then it was
classified as a new well zone.

In the objective functions, maximum pumping rates in the well zones and
maximum groundwater inflow into the river were considered. Conceptual
models in GMS-MODFLOW (Groundwater Modelling System, Aquaveo) were
developed to simulate the river-aquifer exchanges for the lower part of
the River Ain, France and these were coupled with optimization models
available in the MATLAB platform. To implement the coupling between the
two systems, MATLAB scripts were developed which carried out the
following tasks, (i) run the optimization model with Modflow as cost
function, (ii) read and write the decision variables from the `.h5'
input file for Modflow, (iii) run the simulation model, extract and
post-process the output of simulation model to calculate a cost for
optimization model.

Additionally, the sets of optimized costs (Pareto fronts) were compared
using various criteria (convergence, diversity, and uniformity) and
metrics (hypervolume, spread, and IGD). Decision variables corresponding
to these Pareto fronts were also compared in terms of spatial and
quantitative distribution in different model domains.

\textbf{3.1 Groundwater Model Development} :

The groundwater flow modeling was performed using MODFLOW. GMS 10.4
based on conceptual model approach was adopted to create different
geospatial data-based input layers for defining the surface recharge,
boundary conditions, discharge wells, observation wells, and hydraulic
conductivity of the area and to build the groundwater model. The
piezometric surface of groundwater was created using hydrograph data of
280 wells in the study area.

The topography, top, and bottom of the study area, was found in the
range of 240 m to 550 m where SRTM data was used to get the surface
elevation of the area. The bottom surface was prepared with the help of
well log data obtained by BRGM. The two-layered model was developed for
the underlying sediments where each layer was assumed to be horizontal,
homogeneous, and isotropic. The mean thickness of the layer was taken as
25 m. The initial values of hydraulic conductivity were taken from
0.0018 m/s for the older sediment and 0.003 m/s for the younger
sediment. The horizontal hydraulic conductivity was estimated from the
pumping well test data and the literature data obtained by BRGM. Thus,
the obtained data have been used as initial distributions that were
subsequently modified during calibration of the numerical model to
achieve the best fit between the simulated and the observed data.
Specific yield (Ss) values for the alluvial deposits were found in the
range from 1 to 17\%.

\emph{Boundary Conditions:} The model was developed by defining two
types of external model boundaries i.e. constant head and constant flow
boundaries. Figure-3 shows the different types of domains that were
considered to perform the groundwater modeling. The figure shows that in
the Domain-1, eastern \& western sides of the model domain were defined
based on watershed divide line i.e. no-flow boundary respectively,
whereas in Domain-2 alluvium plain of the Ain river were considered and
the eastern \& western side of the model domain was defined based on
constant flow boundaries. The constant flow values were referred to from
the published report of BRGM (Ref.). Whereas, the model Domain-3 was
chosen the same as Domain-1 with one modification in the lower eastern
part. The Rhone river was introduced in the model domain through the
constant head boundary. This modification helped to find the impact of
incorporating the Rhone river in the modeling area.

\emph{Piezometer and River water level data:} In this study, a total of
15 piezometers were available, where a few piezometers had data for the
period of 2008 to 2010 and the remaining piezometers covered the data
from 2002 to 2015. The data thus collected, showed that the piezometer
doesn't show very high fluctuation in values. This variation ranged from
3.7m to 1.2m, which showed that the groundwater table, and
correspondingly groundwater scenario, is very stable in the region.
Further, the river water level data, which was measured at 5 locations
of the study area, was obtained from Banque Hydro France. The data from
2002 to 2015 was also examined to understand the trend of data. This
analysis showed a fluctuation range of 2.54 m at Chezy and 1.67m at
point de'Ain on the river Ain. Data of other tributaries such as
Albarine also showed a maximum fluctuation of 1.36m.

\emph{Recharge \& Evapotranspiration:} Rainfall was considered as the
source of groundwater recharge. Rainfall data from the Météo France
database was used to calculate the recharge input values and applied
uniformly over the polygons constructed in the model domain. The
different recharge polygons were developed based on a land-use map
created from satellite imageries (Figure 6). Six categories were taken
for classification purposes such as water body, agriculture, fallow
land, built-up, forest/vegetation, sand. Supervised classification was
performed to create the land-use map where the fusion of Satellite
imageries of Sentinel 1 and Landsat 8 was used. In Sentinel-1 bands, VV
and VH were used whereas in Landsat 8 bands 1 to 7 were used for
classification purposes. Random Forest model was used as classifier and
training accuracy, test accuracy, and kappa coefficient were 99.98\%
98.45\%, and 97.77\% respectively. The initial value of recharge was
taken as 10\%, 50\%, 50\%, 80\%, and 60\% for built-up, agriculture,
vegetation, sand, and fellow land of total precipitation i.e. 1650
mm/year. Estimated evapotranspiration was taken as 638 mm/year with an
extinction depth of 2 m. The potential evapotranspiration was considered
uniformly distributed over the study area. For the future forecast, 10
years' monthly average of rainfall were taken and corresponding recharge
values were used.

\emph{Water Demand and Supply:} Water consumption in the study area is
mainly done by agriculture, the rest is used for domestic purposes and
some for industrial usage. The discharge details about the wells were
computed on the basis of agricultural and domestic demands of the area.
The field survey was carried out to collect the information of water
demand as per the cropping pattern and domestic water consumption. The
total agriculture water demand was identified through the classification
of satellite imagery data for three seasons. The groundwater abstraction
varies seasonally as well as yearly due to variable demands for
irrigation. The quantitative data, available from 2002, are obtained
from three sources - the Rhône-Méditerranée-Corse Water Agency (AERMC),
the Directorate Department of Agriculture and Forestry (DDAF), and the
Association Syndicale of Ain Irrigation (ASIA). The groundwater
abstraction is divided into three uses (Table 1). Thus, the total
catchments identified are 372, representing a total annual volume of
40,130,104 m\textsuperscript{3}. which varies up to 27
Mm\textsuperscript{3} in some conditions. The total water requirement
computed by this method is found to be 106.78 mcm, which is 5\% more
than the value available by administrative authority i.e., 101.43 mcm.

Calibration and Validation: A regional groundwater flow model was
constructed and calibrated to the transient-state condition with a
stress period of four months. In the calibration of the model, the value
of recharge and boundary inflow was taken. Calibration of the model was
performed based on computed and observed values of groundwater head at
20 evenly distributed points in the study area. Based on the
availability of data, the model was calibrated from 2008 to 2010 based
on all piezometers and further from 2010 to 2012 on the basis of
remaining wells. The further model was validated based on data from 2012
to 2015. In the calibration processes, the groundwater head values which
were computed by the model and observed at the observation wells were
analyzed at 95\% confidence level (Figure 4) at four different
locations. The mean square error was computed for the calibration
process. The differences between observed and computed values were found
less in the middle part of the valley whereas on the boundary of the
model domain they were not found within the range of the 95\% confidence
level due to inaccuracy for defining the head flux boundaries.
Therefore, head flux boundary values were modified within the range of
10\%, and the effect on the error was observed. Initially, the constant
flow boundary was calibrated in the study state condition to incorporate
the groundwater in-flow from the adjacent aquifer. The recharge rate was
calibrated with the help of PEST (a popular parameter estimation
program) in the MODFLOW package. Final outputs of the model were
obtained as water table contour maps and directions of groundwater
movement together with mass water balances for the model domain, and
river-aquifer exchanges. Further, river-aquifer exchanges were
calculated with the help of a calibrated model.

\textbf{3.2. Optimization Problem}

The objective of the optimization model is to maximize the water
withdrawal rate from the aquifer (Total discharge, Q) and water gain of
the river from the aquifer (Leakage out, L) while having the least
drawdown at the wells in the domain. The leakage out represents the rate
of discharge of water from the aquifer to the river. The pumping rate is
assumed constant over well zones and specified time steps. The water
withdrawal rate from the well zones, leakage rate out of the aquifer
into the river, and drawdown at the wells, of the last time step, are
considered for optimization. The domain contains different zones of the
wells, based on the distance from the river and municipal zone data. To
achieve the objective, the discharge in these zones are adjusted by the
following optimization algorithm:

\begin{equation}
\text{Maximise}\left\{\begin{matrix}\sum_{i\in R}L_{i}-P\\
\sum_{i=1}^{n_{z}}{N_{i}*Q_{i}}-P\\
\end{matrix}\right.\ \nonumber \\
 \end{equation}\begin{equation}
\left(Q_{i}\right)_{\text{lb}}\leq Q_{i}\leq\left(Q_{i}\right)_{\text{ub}}\nonumber \\
 \end{equation}

where, L\textsubscript{i} = rate of leakage out of the aquifer to the
river; Q\textsubscript{i} = rate of discharge of
i\textsuperscript{th}well zone; R = set of the river grid cells with
leakage out; n\textsubscript{z} = total number of well zones;
N\textsubscript{i} = number of wells in the i\textsuperscript{th} zone;
P = penalty imposed due to drawdown constraint violation;
(Q\textsubscript{i})\textsubscript{lb} = lower bound of the discharge
for the i\textsuperscript{th} well zone;
(Q\textsubscript{i})\textsubscript{ub} = upper bound of the discharge
for the i\textsuperscript{th} well zone.

\begin{equation}
\mathbb{d}_{\text{dist}}=\sqrt{\sum_{i\in W}\left(d_{i}-d_{\text{threshold}}\right)^{2}}\nonumber \\
 \end{equation}\begin{equation}
P=C_{\text{model}}*d_{\text{dist}}\nonumber \\
 \end{equation}

Where, d\textsubscript{i} = drawdown at i\textsuperscript{th} well; W =
set of wells; d\textsubscript{threshold} = threshold value of drawdown,
taken as 2 m; d\textsubscript{dist} = distance-based metric for
drawdown; C\textsubscript{model} = constant to amplify the penalty. The
value of C\textsubscript{model} is calculated such that the penalty is
of the same order as that of the cost. Both the cost has values of the
order 10\textsuperscript{5} and from the multiple (500) random model
evaluations, an expected value of d\textsubscript{dist} is obtained.
Consequently, the value of the C\textsubscript{mdoel} is obtained as
10\textsuperscript{5}/ E(d\textsubscript{dist}).

\textbf{3.3. Optimization Techniques:}

Four techniques viz. MoGA, MOPSO, PS, and MOEA/D were compared on three
groundwater model domain sizes. MATLAB provides inbuilt functions for
the implementation of MoGA and PS. The following subsections discuss the
details of the optimization algorithms.

\textbf{3.3.1 Multi-objective Genetic Algorithm:} In MATLAB,
`gamultiobj' is the function used for the Multi-objective GA. This
function uses a controlled and elitist genetic algorithm (The MathWork
Inc.), which is a variant of NSGA-II (D. Kalyanmoy, 2001). This
algorithm increases the diversity by favoring a variety of individuals
even if they have a lower fitness value. After the initialization,
generic operations (crossover, mutation, and selection) are iteratively
performed. Consequently, crossover fraction, mutation rate, number of
generations, and population size are major parameters affecting the
performance of the algorithm. For the ANN model, these parameters were
set to 0.85, 0.015, 2000, and 200 respectively, and for the simulation
model, these were set as 0.85, 0.015, 150, and 25 respectively.

\textbf{3.3.2 Multi-objective Particle Swarm Optimization:} This
heuristic algorithm is inspired by the movement of a bird flock. Each
bird or particle has a position and local velocity in the feasible
solution domain. The exploration of the domain by these particles is
governed by both local and global velocity, which are computed based on
the personal best and the global best solution obtained in the domain.
In Coello et. al. (2004), the most widely used implementation of MOPSO
introduced a mutation operator that enriches the exploration capability
of the algorithm. V\selectlanguage{ngerman}íctor Martínez-Cagigal (2020) provides a MATLAB
implementation of Coello et. al. (2004), which is used in this work. The
c1 and c2 were set to 2 for both ANN and simulation models. Maximum
iterations and number of particles, for the ANN model, were set to 1000
and 100 respectively. For the simulation model, these were 150 and 25.

\textbf{3.3.3 Pareto Search:} The Pareto search algorithm by MATLAB,
finds the non-dominated solution by the use of Pareto search in a set of
points (archives and iterates). The algorithm uses the poll to find
better solutions and if no better solution is available, then in the
next iteration half multiply the mesh size. Theoretically, the algorithm
converges to points near the true Pareto front (The MathWork Inc.). In
MATLAB, \emph{paretosearch} function is used for its implementation.
Also, it takes fewer function evaluations than \emph{gamultiobj} and it
performs better or equal, in comparison to NSGA-II if there are no
non-linear constraints (The MathWork Inc.). No parameter tuning is
required for the Pareto search algorithm, though stopping criteria based
upon tolerance and time can be adjusted.

\textbf{3.3.4 Multi-objective Evolutionary Algorithm Based on
Decomposition:} MOEA/D solves a multi-objective by decomposing it into
multiple scalar sub-problems and solving them simultaneously. The
solution of these sub-problems is evaluated based on its neighboring
sub-problems, which makes its computational cost lower in each
generation in comparison to NSGA-II (Zhang et.al. 2007). Polynomial
mutation (order n) is used as the mutation function. Several
sub-problems, maximum iterations, percentage of the neighborhood, and
mutation rate are the major parameters for MOEA/D. For the final ANN
model, these parameters were set to 1500, 100, 40\%, and 1/n
respectively. For the simulation model, these were set to 100, 20, 40\%,
and 1/n respectively.

The use of these algorithms is computationally expensive with the
groundwater simulation models, which require few seconds for a single
run and require days to complete the single solution. Therefore, an ANN
model is also developed to speed up the evaluation of cost function. The
efficiency ANN-Optimization model is majored with the performance of the
actual simulation-optimization algorithm by performing a feasible number
of simulation runs. The ANN model was trained with a dataset generated
by simulation models itself. The parameters of optimization algorithms
are set to improve the solution set at the cost of increasing the number
of simulations which is accommodated by the ANN model.

\textbf{3.4. ANN-Optimization model}

ANN is a widely used technique to map non-linear and black box functions
for fast evaluations. Feed Forward Neural Network (FFNN) is the most
commonly used ANN. It consists of multiple parallel layers of memory
units (neurons). Each layer of neuron is fully connected with its
adjacent layer and the strength of the connection is defined as its
weight. The backpropagation technique is used to find these weights such
that error between the actual and predicted values, is minimized.

For this work, nearly 10\textsuperscript{4} random data points, with
zone-wise discharge as input and leakage out (L with penalty) and total
discharge (Q with penalty) as output, were generated using a calibrated
groundwater simulation model. The total data is divided into three
subsets: training (70\%), validation (15\%), and testing (15\%). An FFNN
with a single hidden layer and enough hidden layer size (neurons) is
capable of any input-output mapping (The MathWork Inc.). But the ANN
training time also increases with an increase in the size of the
training dataset and layer size (neurons). This makes the combination of
large training data set and layer size infeasible. Different layer sizes
(20, 10, 25) were tested and the best of the ANN model was selected
based on R-square and Mean Absolute Percentage Error (MAPE) of Leakage
out because the total discharge predictions were accurate. Along with
these matrices, the values of Maximum Absolute Error (MXAE) and Root
Mean Squared Error (RMSE) are shown in Table 3.

As is evident from the table, the ``Model 2'', under leakage out, is
having a good R-square value (near to 0.97) and a high MAPE value (near
to 25\%). This means that even though the actual and predicted output
are correlated, there is a significant deviation existing from the
actual model. Even if the optimum points for Model 2 (ANN) are derived
using evolutionary algorithms, the correct value of decision variables
could not be obtained. It makes the use of this ANN model unsuitable for
the optimization model for Model 2. Thus, final optimization using ANN
is performed for ``Model 1'' and ``Model 3'' only.

\textbf{3.5. Performance Metrics:}

In multi-objective optimization algorithms, the solution set generated
approximates the actual Pareto front. The approximated fronts are
compared based on convergence, diversity, and uniformity. Various
performance matrices are available that account for one, two, or all of
the three aspects of the Pareto fronts. These metrics can be unary or
binary. In this work, six metrics are used, namely, epsilon, spread,
generalized spread, generational distance, inverted generational
distance, and hypervolume (Table 4). The matrices require a reference
point or reference Pareto front. The non-dominated set from the union of
the points obtained from the compared Pareto front is used as the
reference Pareto front and a value lower than that of the anti-utopia
point of this reference front is taken as a reference point (1e4, 5e4).
A detailed comparison of various performance metrics can be found in
Zitzler et. al. (2003).

The unary e-indicator metric represents the smallest distance that an
approximate Pareto front must be translated to completely dominate the
reference Pareto set (Kollat et. al., 2005). The average Euclidean
distance between the reference Pareto set and the approximate Pareto
solutions is called the Generational distance (GD). Both, epsilon and
generational distance measure the convergence. The lower the value of
epsilon and generational distance, the better is the convergence. The
diversity of the Pareto front can be compared based on the distribution
of the solution set (uniformity) and its extent. Spread and generalized
spread quantifies the non-uniformity of approximate Pareto front. Small
the value of these matrices indicates, a better and diverse set of
approximate Pareto front. The inverted GD is the average distance
between each member of the reference Pareto front and the approximate
Pareto front. It is highly influenced by the distribution of the
approximate Pareto front. Hypervolume measures the size of the space
enclosed by the approximate Pareto front. Both, inverted generational
distance and hypervolume are a unary metric that considers convergence
and diversity. In addition, hypervolume is the most widely used
performance metric (Riquelme et. al. 2015). Auger A. et. al. (2009) have
provided a method to choose the value of reference point for the
calculation of hypervolume. Unlike inverted generational distance, the
higher the hypervolume better is the combined effect of convergence and
diversity of the approximate Pareto front.

\textbf{4. Results:}

\textbf{4.1 Simulation-optimization model}

Each combination of the simulation model and optimization technique was
evaluated up to 5000 points depending upon the optimization algorithm
and its parameters. In all the three models, the Pareto fronts by MOPSO
were more converged, whereas the results obtained from GA were more
diverse for model 1 and model 3 (Fig. 5). The Pareto search and MOEAD
were more converged in comparison to GA, but they lack diversity. Only
in model 2, the diversity of the PF of MOPSO outperforms that of GA. The
solutions of GA are highly diverse in model 1 and model 3, even at a low
number of simulations runs. Other than the optimization algorithm, the
Pareto fronts trends are determined by the model demarcation and its
boundary conditions.

To understand the impact of different model domains on optimal outputs,
the best-performed MOPSO front was studied. Pareto fronts show that
Domain-1 is giving a high value of leakage out in comparison of the
optimal value of groundwater discharge through wells whereas the output
of Domain-2 is higher than Domain-3 ex. for discharge value of 2000000
m\textsuperscript{3}/day leakage out in Domain-1, Domain-2 and Domain-3
are as follows 153560 m\textsuperscript{3}/day,134940
m\textsuperscript{3}/day and 100800 m\textsuperscript{3}/day. The Pareto
fronts (Fig. 5) also depict the effect of a model domain on the optimal
discharge and river gain relation. For instance, in the case of MOPSO,
an increase of discharge by 50,000 m\textsuperscript{3}/day (from
200,000 m\textsuperscript{3}/day to 250,000 m\textsuperscript{3}/day)
leads to a decrease in river gain by 15,820 m\textsuperscript{3}/day,
37,657 m\textsuperscript{3}/day, and 19,535 m\textsuperscript{3}/day for
the model domain 1, 2, and 3 respectively. Correspondingly, the Pareto
search demonstrates a decrease of river gain by 22,278.9
m\textsuperscript{3}/day, 40,116 m\textsuperscript{3}/day, and 20,457.4
m\textsuperscript{3}/day for the previously mentioned discharge
increase.

All the six metrics, for the three model domains and four optimization
techniques, are shown in Figure 6. The combined effect of diversity and
convergence from inverted GD and hypervolume also suggests the
superiority of MOPSO, in all three domains, over other optimization
techniques. Further, the epsilon metric and generational distance
suggest that convergence of MOPSO is best. Both the parameters, epsilon,
and generational distance, show conflict in the comparison of other
optimization techniques. In terms of uniformity of solutions, spread and
generalized spread suggests the superiority of GA for model 1. In model
3, even though the Pareto front of GA is diverse, the solutions are not
uniform due to which the spread of GA is high. For model 3, MOPSO is
much more uniformly distributed and hence has a low value of the spread.

\textbf{4.2 ANN-Optimization model}

MOPSO showed significant improvement in terms of convergence, diversity,
and uniformity of the solutions (Figure 7). MOEAD showed improved
diversity and uniformity whereas convergence was still similar to the
Pareto search (without ANN). Although GA showed no major improvement in
the convergence, the uniformity of Pareto solution distribution is
increased. The Pareto search solution without the use of ANN can be
taken as the reference for the comparison (Fig. 5 and Fig. 7). The
diversity and uniformity of the GA and MOPSO are close and comparable.

Spread and Generalized spread depict that the MOPSO and GA have a
uniform distribution of the solutions. Whereas the MOEAD shows a polar
and non-uniform distribution of solutions. The MOPSO is superior to
other algorithms, in terms of convergence, suggested by epsilon and
generational distance. The combined effect of convergence and diversity,
measured with hypervolume and inverted GD, is also best for MOPSO. GA,
PS, and MOEAD show a close value of hypervolume and inverted GD
indicating a nearly equivalent performance in terms of convergence and
diversity. The percentage increase in hypervolume and decrease in
inverted GD, both show that MOPSO in model 1and MOEAD in model 3, has
significant improvement (Fig. 9). Whereas GA and MOEAD in model 1 show
the least improvement. However, the GA and MOPSO (in Model 3) showed
conflict in the change in values. Inverted GD is less reliable because
it is strongly influenced by the distribution of the approximate Pareto
front. Figure 10(a) and Figure 10(b) shows the box-whisker plot of the
decision variables from the Pareto solutions obtained from each
optimization technique and two (1 and 3) model domain size. It can be
observed that the distribution of decision variables for the solutions
is dependent on the optimization algorithm rather than the model domain
size. Further, in GA the solution decision variables have a smaller
number of outliers compared based on well zone distribution. A similar
trend is shown in MOPSO and Pareto searches, with an increased number of
outliers. MOEAD shows a broader distribution with the length of the box
(25 and 75 percentiles bound) greater than other optimization
techniques. MOPSO solutions have a large number of outliers depicting a
large number of solutions focused in the quartile range. The mean
discharge for GA solutions is less than or close to 1000 for nearly all
of the zones. In Pareto Search, the mean values are close to the bounds
of the box (25 or 75 percentile) of the solutions and the mean discharge
largely varies with the decision variable number (i.e., location of well
zone).

Furthermore, figure 10 provides the spatial variation of the optimized
well discharge. In total, there are 30 well zones in both the domains
and half of them are classified as ``far away'' or ``close'' depending
upon their shortest distance from the river. The average optimal
discharge of faraway wells is observed to be slightly higher than those
close to the river. In Model 3, the value means optimal discharge of
faraway well is 860.8 m\textsuperscript{3}/day and that of wells close
to the river is 832.9 m\textsuperscript{3}/day. However, in Model 1 the
gap between the two values narrows down further, with far away well avg.
discharge as 905.4 m\textsuperscript{3}/day and close well avg.
discharge as 901.2 m\textsuperscript{3}/day.

\textbf{5. Discussions} :

The solutions of GA are highly diverse in model 1 and model 3, even at a
low number of simulations runs. Other than the optimization algorithm,
the Pareto fronts trends are determined by the model demarcation and its
boundary conditions. It shows that adding the Rhone River in the model
alters the groundwater flow correspondingly water budgeting of the area
and gives the least groundwater flow into the river due to excess
pumping. Results show that Model-2 which consists of the shortest area
along with the dominated constant boundary inflow is providing less
groundwater to the river in comparison to Domain-1 which consists of a
larger area in the simulation along with a one-sided watershed boundary.
It can be observed that even though all three domains extract the same
amount of groundwater, the modelled effect on the river varies for each
domain. In addition, the data from the Pareto front suggest that in
model 2, the extraction from the well and the river gain are heavily
related to each other. This relation will not only influence the
simulation results but will also affect the decision-making in
groundwater management. In ANN Optimization model, it is thus evident
that in terms of convergence, the solution points of MOPSO dominate
every other solution, and solution points of GA are dominated by most of
the other points. PS and MOEAD solutions lie in between in the GA and
MOPSO front for both models. MOPSO solutions largely consist of outliers
with either too small or too large box bounds. MOEAD solutions are
largely uniform and cover the lower and upper bound of decision
variables. Regardless of the small difference, the overall trend is:
wells close to rivers extract less groundwater than far away. Another
factor, that can contribute to the high discharge of a close well is
spatial heterogeneity in the aquifer properties.

\textbf{5. Conclusions:}

Three different domain sizes of groundwater model for optimization of
withdrawal and gain in the river from the aquifer, by four different
optimization algorithms were discussed. The groundwater model
development, optimization models, ANN model to reduce evaluation time,
performance metrics for Pareto fronts, and the comparison of the
distribution of optimal decision variables for different optimization
techniques, are also discussed. The ANN model for three domains shows
the difference in accuracy, suggesting a change in domain size and
consequently boundary conditions can alter the performance of ANN.
Results show that boundary conditions and domain size influence the
result of simulation-optimization models. Domain-1 was found more
efficient and was capable to give a higher value of R-A exchange
corresponding to the discharge of pumping wells.

The simulation-optimization model result suggests that GA provided a
diverse set of Pareto solutions even for a low number of evaluations but
the solution of other algorithms was dominating, with MOPSO being the
best. And among the three domains, model 2 consistently implied a
greater interrelation between the two costs, i.e. total discharge and
river gain. In the ANN-Optimization model, MOPSO showed significant
improvement in both diversity and convergence. GA solutions didn't show
major improvement in convergence. This suggested the use of GA to obtain
a tentative solution set that can be obtained without the use of the ANN
model. The use of ANN can significantly improve the performance of MOPSO
by allowing it to have a greater number of evaluations, which was
limited in the raw simulation model due to the infeasible computation
time.

The distribution of optimal decision variables was well distributed for
MOEAD and comparatively limited in other optimization techniques. But
the distribution of decision variables didn't represent the quality of
the solution from perspective of performance of optimization algorithm.
In case of decision-making problems, this study can be helpful to have
good distribution of decision variables, even though they provide a
better solution but not the best one. Apart from statistical
distribution of solutions, the spatial distribution of optimized
decision variables (well discharge) provided a lot of information
regarding the physical interpretation and conceptual verification of the
solutions. These solutions followed an expected trend of higher
discharge at faraway wells and lower discharge at the close wells. From
this study, the influence of domain demarcation, boundary conditions,
and optimization algorithm has been observed on the Pareto fronts and
optimized decision variables. The conceptual verification of the Pareto
solutions was also explained. Besides this, the use of ANN and handling
the constraints of the optimization problem is also considered for a
real-world groundwater problem. Still, several challenges in R-A
exchange-related groundwater decision-making exist. Accurate R-A
exchange modeling, precise surrogate models to reduce S-O time, and
dealing with high heterogeneity of the domain, to name a few.

\textbf{References}

\begin{itemize}
\tightlist
\item
  Dangar, Swarup \& Asoka, Akarsh \& Mishra, Vimal. (2021). Causes and
  implications of groundwater depletion in India: A review. Journal of
  Hydrology. 596. 10.1016/j.jhydrol.2021.126103.
\item
  Anderson MP. (2005). Heat as a Ground Water Tracer. Ground Water,
  43(6), 951-968.
\item
  Asadzadeh, Masoud, Bryan A. Tolson, and Donald H. Burn (2014). ''A new
  selection metric for multiobjective hydrologic model calibration.''
  Water Resources Research 50, no. 9, 7082-7099.
\item
  Atwell BH, MacDonald RB, Bartolucci LA. (1971). Thermal mapping of
  streams from airborne radiometric scanning. Water Resources Bulletin,
  7(2), 228-243.
\item
  Auger, A., Bader, J., Brockhoff, D. and Zitzler, E., (2009). Theory of
  the hyper volume indicator: optimal \selectlanguage{greek}μ-\selectlanguage{english}distributions and the choice of
  the reference point. In Proceedings of the tenth ACM SIGEVO Workshop
  on Foundations of genetic algorithms, 87-102.
\item
  Avent, Brendan, Javier Gonz\selectlanguage{ngerman}ález, Tom Diethe, Andrei Paleyes, and Borja
  Balle (2020). ''Automatic Discovery of Privacy--Utility Pareto
  Fronts.''~Proceedings on Privacy Enhancing Technologies~2020, no. 4,
  5-23.
\item
  Bassi, Mohamed, Eduardo Souza de Cursi, Emmanuel Pagnacco, and Rachid
  Ellaia (2018). ''Statistics of the Pareto front in Multi-objective
  Optimization under Uncertainties.'' Latin American Journal of Solids
  and Structures 15, no. 11.
\item
  Becker MW, Georgian T, Ambrose H, Siniscalchi J, Fredrick K. (2004).
  Estimating flow and flux of groundwater discharge using water
  temperature and velocity, Journal of Hydrology, 296(1-4), 221-233.
\item
  Belakaria, Syrine, and Aryan Deshwal (2019). Max-value entropy search
  for multi-objective bayesian optimization. In International Conference
  on Neural Information Processing Systems (NeurIPS).
\item
  Binois, M., Ginsbourger, D., \& Roustant, O. (2015). Quantifying
  uncertainty on Pareto fronts with Gaussian process conditional
  simulations. European journal of operational research, 243(2),
  386-394.
\item
  Burkholder BK, Grant GE, Haggerty R, Khangaonkar T, Wamper PJ. (2008).
  Influence of hyporheic flow and geomorphology on temperature of a
  large, gravel-bed River, Clackamas River, Oregon, USA. Hydrological
  Processes, 22, 941-953.
\item
  Calandra, Roberto, Jan Peters, and M. P. Deisenrothy (2014). Pareto
  front modeling for sensitivity analysis in multi-objective bayesian
  optimization. In NIPS Workshop on Bayesian Optimization, vol. 5.
\item
  Cao, P., Fan, Z., Gao, R., \& Tang, J. (2017). A
  manufacturing-oriented single point search hyper-heuristic scheme for
  multi-objective optimization. In ASME 2017 International Design
  Engineering Technical Conferences and Computers and Information in
  Engineering Conference. American Society of Mechanical Engineers
  Digital Collection.
\item
  Coello, Carlos A. Coello, Gregorio Toscano Pulido, and M. Salazar
  Lechuga (2004). ''Handling multiple objectives with particle swarm
  optimization.'' IEEE Transactions on evolutionary computation 8.3,
  256-279
\item
  Constantz J. (1998). Interaction between stream temperature,
  streamflow, and groundwater exchanges in alpine streams. Water
  Resources Research, 34, 1609-1615.
\item
  Cristea NC, Burges SJ. (2009). Use of Thermal Infrared Imagery to
  Complement Monitoring and Modeling of Spatial Stream Temperature.
  Journal of Hydrologic Engineering, 14(10), 1080-1090.
\item
  Deb, Kalyanmoy (2001). Multi-Objective Optimization using Evolutionary
  Algorithms, John Wiley \& Sons, Ltd, Chichester, England.
\item
  Emmerich, Michael TM, and André H. Deutz (2018). A tutorial on
  multiobjective optimization: fundamentals and evolutionary methods.
  Natural Computing 17, no. 3, 585-609.
\item
  Graillot D, Paran F, Bornette G, Marmonier P, Piscart C, Cadilhac L.
  (2014). Coupling groundwater modeling and biological indicators for
  identifying river/aquifer exchanges. SpringlerPlus, DOI:
  10.1186/2193-1801-3-68
\item
  Handcock RN, Gillespie AR, Cherkauer KA, Kay JE, Burges SJ, Kampf SK.
  (2006). Accuracy and uncertainty of thermal infrared remote sensing of
  stream temperatures at multiple spatial scales. Remote Sensing of
  Environment, 100, 427-440.
\item
  Herbert C, Caissie D, Satish MG, El-Jabi N. (2011). Study of stream
  temperature dynamics and corresponding heat fluxes within Miramichi
  River catchments (New Brunswick, Canada). Hydrological Processes, 25,
  2439-2455.
\item
  Hernández-Lobato, Daniel, Jose Hernandez-Lobato, Amar Shah, and Ryan
  Adams (2016). ''Predictive entropy search for multi-objective bayesian
  optimization.'' In International Conference on Machine Learning, pp.
  1492-1501. PMLR.
\item
  Höllermann, Britta, and Mariele Evers (2019). Coping with uncertainty
  in water management: Qualitative system analysis as a vehicle to
  visualize the plurality of practitioners' uncertainty handling
  routines. Journal of environmental management, 235, 213-223.
\item
  Horn, D., Demircio\selectlanguage{polish}ğ\selectlanguage{english}lu, A., Bischl, B. et al. (2018) A comparative
  study on large-scale kernelized support vector machines. Adv Data Anal
  Classif 12, 867--883.
\item
  Kalbus E, Reinstorf F, Schirmer M. (2006). Estimating flow and flux of
  groundwater discharge using water temperature and velocity. Hydrology
  and Earth System Sciences, 10, 873-887.
\item
  Kay JE, Kampf SK, Handcock RN, Cherkauer KA, Gillespie AR, Burges SJ.
  (2005). Accuracy of lake and stream temperatures estimated from
  thermal infrared images. Journal of the American Water Resources
  Association, 41, 1161-1175.
\item
  Keery J, Binley A, Crook N, Smith JWN. (2007). Temporal and spatial
  variability of groundwater-surface water fluxes: Development and
  application of an analytical method using temperature time series.
  Journal of Hydrology, 336, 1-16.
\item
  Kollat, Joshua B., and Patrick M. Reed (2005). The value of online
  adaptive search: a performance comparison of NSGAII, \selectlanguage{greek}ε-\selectlanguage{english}NSGAII, and
  \selectlanguage{greek}ε\selectlanguage{english}MOEA. International Conference on Evolutionary Multi-Criterion
  Optimization. Springer, Berlin, Heidelberg.
\item
  Loheide SP, Gorelick SM. (2006). Quantifying stream -- Aquifer
  Interactions through the Analysis of Remotely Sensed Thermographic
  Profiles and In Situ Temperature Histories. Environmental Science and
  Technology, 40(10), 3336-3341.
\item
  Lowry CS, Walker JF, Hunt RJ, Anderson MP. (2007). Identifying spatial
  variability of groundwater discharge in a wetland stream using a
  distributed temperature sensor. Water Resources Research, 43, W10408.
\item
  Marjit R. and Hopfe C.J. (2009). Multi-objective robust optimization
  algorithms for improving energy consumption and thermal comfort of
  buildings.
\item
  MATLAB. (20\_\_).~version 7.10.0 (R20\_\_a). Natick, Massachusetts:
  The MathWorks Inc.
\item
  McDonald MG, Harbaugh AW. (1988). A modular three-dimensional
  finite-difference ground-water flow model, 6. U.S. Geological Survey
  Techniques of Water-Resources Investigations.
\item
  Paran F, Arthaud F, Bornette G, Graillot D, Lalot E, Marmonier P,
  Novel M, Piscart C. (2012). Characterization of exchanges between
  Rh\selectlanguage{ngerman}ône River and groundwater., Research report ZABR Phase 4, AERMC
\item
  Riquelme, Nery, Christian Von Lücken, and Benjamin Baran.(2015)
  Performance metrics in multi-objective optimization. Latin American
  Computing Conference (CLEI). IEEE.
\item
  Salomon, Ludovic, Charles Audet, Jean Bigeon, and Sébastien Le Digabel
  (2018). Review of the quality of approximated Pareto fronts in
  multiobjective optimization, Journées de l'optimisation. In JOpt.
\item
  Sophocleous, M. (2002) Interactions between groundwater and surface
  water: the state of the science,~Hydrogeology
  Journal,~10\textbf{,~}52--67.
\item
  Víctor Martínez-Cagigal (2020). Multi-Objective Particle Swarm
  Optimization (MOPSO) (The MathWorks Inc.), MATLAB Central File
  Exchange. Retrieved October 1, 2020.
\item
  Wawrzyniak V, Piégay H, Poirel A. (2012). Longitudinal and temporal
  thermal patterns of the French Rhône River using Landsat ETM+ thermal
  infrared images. Aquatic Sciences, 74(3), 405-414.
\item
  Wei Mao, Jinzhong Yang, Yan Zhu, Ming Ye, Jingwei Wu (2017) Loosely
  coupled SaltMod for simulating groundwater and salt dynamics under
  well-canal conjunctive irrigation in semi-arid areas, Agricultural
  Water Management, 192, 209-220.
\item
  Westhoff MC, Savenije HHG, Luxemburg WMJ, Stelling GS, van de Giesen
  NC, Selker JS, Pfister L, Uhlenbrook S. (2007). A distributed stream
  temperature model using high resolution temperature observations.
  Hydrology and Earth System Sciences, 11, 1469-1480.
\item
  Zhang, Qingfu, and Hui Li (2007). MOEA/D: A multiobjective
  evolutionary algorithm based on decomposition. IEEE Transactions on
  evolutionary computation 11.6, 712-731
\item
  Zitzler, Eckart, Lothar Thiele, Marco Laumanns, Carlos M. Fonseca, and
  Viviane Grunert Da Fonseca. Performance assessment of multiobjective
  optimizers: An analysis and review. IEEE Transactions on evolutionary
  computation 7, no. 2 (2003): 117-132.
\end{itemize}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/study-Area/study-Area}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/methodolody/methodolody}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/03-domains/03-domains}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/04-calibration-graph/04-calibration-graph}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/05-pareto-front-comparison/05-pareto-front-comparison}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/06-metrices-for-comparison-of-optim/06-metrices-for-comparison-of-optim}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/07-pareto-front-comparison/07-pareto-front-comparison}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/08-metrices-of-comparison/08-metrices-of-comparison}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/09-percentage-increse-hypervlom/09-percentage-increse-hypervlom}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/10/10}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/11/11}

\end{center}
\end{figure}

\textbf{Hosted file}

\verb`Tables.docx`  available at \url{https://authorea.com/users/435384/articles/538448-optimization-of-groundwater-pumping-and-river-aquifer-exchanges-for-management-of-water-resources}

\selectlanguage{english}
\FloatBarrier
\end{document}

