Groundwater monitoring plays a significant role in groundwater management. This study presents an optimization method for designing groundwater-level monitoring networks. The proposed design method was used in the Eshtehard aquifer, in central Iran. Three scenarios were considered to optimize the locations of the observation wells: (1) designing new monitoring networks, (2) redesigning existing monitoring networks, and (3) expanding existing monitoring networks. The kriging method was utilized to determine groundwater levels at non-monitoring locations for preparing the design data base. The optimization of the groundwater monitoring network had the objectives of (1) minimizing the root mean square error and (2) minimizing the number of wells. The non-dominated sorting genetic algorithm (NSGA-II) was applied to optimize the network. Inverse distance weighting interpolation was used in NSGA-II to estimate the groundwater levels while optimizing network design. Results of the study indicate that the proposed method successfully optimizes the design of groundwater monitoring networks that achieve accuracy and cost-effectiveness.

Groundwater monitoring plays an important role in collecting data to assess changes in environmental processes in groundwater resources contamination. Loáiciga et al. (1992) classified groundwater monitoring design into hydrogeological approaches and statistical approaches. Hydrogeological approaches use hydrogeologic information and expert judgment to assess and design groundwater monitoring networks. Statistical approaches acknowledge the uncertainty associated with human knowledge of the underlying hydrogeology and treat aquifer properties as random or spatially correlated variables. Statistical approaches include simulation-based approaches, variance-based approaches, and probability-based approaches.

McKinney & Loucks (1992) developed a network design algorithm to improve the reliability of groundwater simulation model predictions. Hudak (2006) reported a Monte Carlo (MC) physics-based simulation approach to locate detection wells in aquifers beneath landfills. Yang et al. (2008) used the kriging standard deviation as a criterion for determining the density of groundwater-level monitoring networks. Dokou & Pinder (2009) created an optimal search strategy that identified contamination sources by the least number of water quality samples. The strategy included a MC stochastic groundwater flow and transport model, a predetermined set of potential source locations, and a Kalman filter which updated the simulated contaminant concentration field using contaminant concentration data. Mogheir et al. (2003) introduced the spatial structure by means of a trans-information and correlation model for groundwater quality variables. The trans-information model was superior to the correlation model in describing the spatial variability (structure) of groundwater quality variables. Wu et al. (2006) compared the MC simple genetic algorithm (MCSGA) and noisy genetic algorithm (NGA) in the design of cost-effective sampling networks considering uncertainties in the hydraulic conductivity. Both methods combined the genetic algorithm (GA) with a numerical flow and transport simulator and a global plume estimator to optimize the sampling network for contaminant plume monitoring. Dhar & Patil (2012) designed groundwater quality monitoring networks under epistemic uncertainty considering spatiotemporal pollutant concentrations as fuzzy numbers. The proposed methodology incorporated fuzzy ordinary kriging (FOK) within the decision model formulation for spatial estimation of contaminant concentration values. The non-dominated sorting genetic algorithm (NSGA-II) was used to solve a design model. Results of the study showed the applicability of the proposed methodology for network design under epistemic uncertainty.

In recent decades, several groundwater monitoring studies have turned to machine learning approaches to add or remove monitoring stations. Data mining is an analytic process to explore data by consistent patterns and/or systematic relation between variables (Fallah-Mehdipour et al. 2013a, 2013b, 2014; Orouji et al. 2013, 2014; Akbari-Alashti et al. 2014; Aboutalebi et al. 2015, 2016a, 2016b, 2016c; Soleimani et al. 2016a, 2016b; Bozorg-Haddad et al. 2017). Asefa et al. (2004) described a methodology based on support vector machines to design monitoring networks. Khader & McKee (2014) applied a regression vector machine (RVM) for groundwater monitoring network design. Their RVM method employed a MC simulation process to capture the uncertainties in recharge, hydraulic conductivity, and nitrate reaction processes. This paper presents an optimization method to design reliable and efficient groundwater monitoring networks. The method has as objectives reducing costs and increasing the groundwater-level monitoring accuracy.

Interpolation methods are used to determine the values of a variable at any point when values of the variable are available at a set of sampling points in a region. Some form of weighted average of the values of a variable (or variables) at surrounding points is applied to calculate the value at the point where the values are unknown. Therefore, the estimation at a location where a spatial variable is unknown is usually given by the weighted average of the nearest neighbors (see e.g., Tobler 2008):
(1)
where = unknown value of at time t, = available measurements of a spatial variable at location j and time t, and = weight applied to , which differs according to the interpolation method used. Spatial interpolation methods belong to two main categories: (1) deterministic (e.g., inverse distance weighting (IDW), splines, radial basis functions, etc.) and (2) geostatistical (e.g., kriging, hierarchical models, copula, etc.). The former methods use a mathematical function to calculate values at unknown locations and provide deterministic estimates. The latter methods provide probabilistic estimates of a variable and its variance of estimation at points where measurements do not exist (Loáiciga et al. 2010).
IDW estimation at an unknown location is obtained with the weighted average of all available measurements. The weights are proportional to the inverse of the distances between the location of interest and those where measurements of the estimated variable are available. The closest available observations are weighted more heavily:
(2)
where = weights applied to the distance between and . The weights usually are chosen as a power function of the Euclidean distance between two spatial points, or , where p = 2 (Bivand et al. 2008).
Kriging is a common method of geostatistics used for spatial interpolation. There are several types of kriging, such as simple, ordinary, universal, block, regression, and cokriging (see Cressie 1991). A brief description of kriging is defined here. Variogram is variance of the difference between values of variable at two locations which is used for spatial interpolation by kriging. The first step of kriging is estimating the empirical semivariogram:
(3)
where = the semivariogram value at distance h, = the total number of the variable pairs separated by a distance h, Z(xi)= the value of the variable at location xi. A parametric or nonparametric model such as Gaussian, spherical, exponential, linear models is fitted to the empirical semivariogram.
In the second step of kriging, the weights in the general interpolation (Equation (1)) are estimated from the following system of equations:
(4)
where μ= Lagrange multiplier and = semivariogram between and .
The minimum squared error of estimation (kriging variance) is a measure for the accuracy of the estimates, given by:
(5)
where = kriging variance at point . Details about kriging can be found, among others, in Cressie (1991).

Multi-objective evolutionary algorithms (MOEAs) were divided into two categories (Deb et al. 2002): (1) non-elitist MOEA and (2) elitist MOEA. Elitism in evolutionary algorithms is the use of superior solutions of a past generation to the next generation solutions. One of the most common methods of elitist MOEAs is the NSGA-II, an improved version of the NSGA (Srinivas & Deb 1994).

The NSGA-II is one of the most common and effective algorithms for solving multi-objective problems. It is a random-based search algorithm and a variant of the GA.

The NSGA-II applies selection, crossover, and mutation operators. The crossover operator recombines the members of a population to make a new population. The mutation operator is applied to manipulate the population's members. The NSGA-II begins with the random generation of a population, which is subsequently sorted based on non-domination into several Pareto fronts. The first Pareto front is a completely non-dominated set whose members are not dominated by the members of other fronts. The members of the second front dominate those from subsequent fronts, but are dominated by the members of the first front, and so on and so forth.

Each member of a front is assigned a fitness value or rank (Deb et al. 2002). For example, individuals from the first front have rank 1, and those from the members of the next front have a rank equal to 2, etc. The crowding distance is a parameter of the NSGA-II that measures the distance of a member to its neighbors and ensures diversity in a population.

As soon as the rank and the crowding distance for all of the members of all fronts are determined, then parents are selected from the population by using binary tournament selection based on the rank and the crowding distance. An individual having smaller value of rank or greater crowding distance is selected (Figure 1). The selected population uses crossover and mutation operators to generate offspring. The current population and current offspring are sorted again and the best individuals are selected according to their rank and crowding distances. Figure 2 shows a flowchart of the NSGA-II.

Figure 1

Process of sorting in the NSGA-II (t = number of generation, P = parent, O = offspring, and Fi = front of rank i).

Figure 1

Process of sorting in the NSGA-II (t = number of generation, P = parent, O = offspring, and Fi = front of rank i).

Close modal
Figure 2

Flowchart of the NSGA-II method.

Figure 2

Flowchart of the NSGA-II method.

Close modal

The study area is the aquifer of the Eshtehard plain, which covers an area equal to 235 km2 in north central Iran (see a map in Figure 3). The plain is surrounded by the Qazvin plain and the Tehran-Karaj plain to the west and east, respectively. The study area features an arid climate, thus groundwater is the main resource of water for residents. The uncontrolled exploitation of groundwater in the plain has produced a declining groundwater table. There is a need for an optimal monitoring network to characterize aquifer conditions as groundwater is mined to meet several water-supply functions. The existing groundwater monitoring network in the region has 18 wells. A four-year period (2009–2013) of recorded data (historical) of existing observation wells was chosen for the study.

Figure 3

Location of the study area.

Figure 3

Location of the study area.

Close modal

First, a groundwater data base was prepared prior to designing a monitoring network and choosing optimal locations of its observation wells within the study area. Next, an optimal groundwater monitoring network was designed.

Data base

Available records of groundwater level were used for estimating the groundwater level over the entire Eshtehard aquifer using kriging.

Optimizing model

The two objective functions of the groundwater monitoring network design method are:

  • 1.

    minimizing the number of observation wells in the area; and

  • 2.
    minimizing the root mean square error (RMSE) between observational and estimated values at all locations within the aquifer:
    (6)
    (7)
    in which f1 = the first objective function (number of wells), f2 = the second objective function (RMSE), = observational groundwater level of point i at the end of period t, = groundwater-level estimation of point i at the end of period of t, T= number of total time periods, N= number of groundwater-level estimation points.

Estimated values of groundwater level were obtained for all the potential sampling locations using IDW in the optimization phase. The RMSE of the network is determined based on Equation (7). The estimated values with IDW and the observational values of the potential points are used to construct the first objective function (minimization of the number of wells). Figure 4 shows a flowchart of the methodology.

Figure 4

Flowchart of groundwater monitoring network design process.

Figure 4

Flowchart of groundwater monitoring network design process.

Close modal

Optimization algorithm

The NSGA-II was selected for solving the multi-objective optimization of the groundwater-level monitoring network. Three scenarios of monitoring network designing were considered:

First scenario

Redesigning of groundwater monitoring network in the study region. The aim of solving this scenario was finding the best location of observation wells over all areas of the aquifer. In this scenario, the groundwater monitoring network was redesigned regardless of the position of available observation wells and the best locations of the wells and RMSE of the network were obtained for each observation well.

Second scenario

Selecting the best set of wells among existing observation wells in the study area. This scenario allows the use of several existing observation wells among the existing monitoring locations to achieve the optimization objectives. For instance, 18 observation wells exist in the research area, of which five of the 18 wells were chosen using the joint RMSE. Therefore, this scenario redesigns an existing monitoring network by keeping active a subset of the existing monitoring wells.

Third scenario

Adding extra wells to the observation wells of the network in the research area. This scenario expands an existing monitoring network considering the joint RMSE.

A data base of groundwater level in the Eshtehard aquifer, central Iran, was established using kriging. Groundwater-level data from 18 observation wells for a four-year period measured monthly was used for interpolation purposes. The NSGA-II algorithm was applied for solving the monitoring network optimization problem. The number of iterations, population size, cross-over, and mutation probabilities were set equal to 1,000, 50, 0.7, and 0.2, respectively, in the NSGA-II. A maximum of 30 wells was designated for the optimized monitoring network. The summary of results is as follows.

A – Scenario 1

The NSGA-II found the optimal groundwater monitoring locations in the aquifer and in the corresponding results were represented as a Pareto front (see Figure 5). Due to the randomized nature of the algorithm's solution, three separate runs were performed to find representative results. The calculated Pareto fronts shown in Figure 5 imply very close results for the three runs, which, in turn, shows the reliable convergence of the NSGA-II over several runs.

Figure 5

Pareto fronts of the monitoring network design problem for three different runs (first scenario).

Figure 5

Pareto fronts of the monitoring network design problem for three different runs (first scenario).

Close modal

Figure 6 shows a sample of 15 optimized groundwater monitoring wells obtained in the first run of the first scenario. It is clear in Figure 6 that the wells have a suitable distribution so that they cover all areas of the aquifer.

Figure 6

Selected wells of the monitoring network (15 wells, first scenario).

Figure 6

Selected wells of the monitoring network (15 wells, first scenario).

Close modal

B – Scenario 2

Under this scenario the monitoring network was redesigned from the existing observation wells in the Eshtehard area. Results of this scenario are presented in Figure 7 in the form of three Pareto fronts obtained in three separate runs. Notice the closeness of the Pareto fronts, which means the NSGA-II algorithm converged to almost the same solution in all runs.

Figure 7

Pareto front of the monitoring network design problem for three different runs (second scenario).

Figure 7

Pareto front of the monitoring network design problem for three different runs (second scenario).

Close modal

A sample of optimized groundwater monitoring locations is presented in Figure 8.

Figure 8

Selected wells of the monitoring network for (15 wells, second scenario).

Figure 8

Selected wells of the monitoring network for (15 wells, second scenario).

Close modal

C – Scenario 3

Under this scenario 18 available observation wells within the aquifer area were kept and the monitoring network was expanded with a few more wells. The results of the optimization are depicted as Pareto fronts in Figure 9, in which the closeness of the three Pareto fronts is evident. Given the randomized algorithmic search it is clear that the similarity of the results of the three runs imply adequate convergence to a near globally optimal solution.

Figure 9

Pareto front of the monitoring network design problem for three different runs (third scenario).

Figure 9

Pareto front of the monitoring network design problem for three different runs (third scenario).

Close modal

A sample of the optimized groundwater monitoring network under Scenario 3 is portrayed in Figure 10. The chosen monitoring wells are suitably distributed throughout the aquifer.

Figure 10

Selected wells of the monitoring network (25 wells, third scenario).

Figure 10

Selected wells of the monitoring network (25 wells, third scenario).

Close modal

Results of the first runs of the three scenarios are presented in Figure 11 as Pareto fronts. The flexible nature of the first scenario, compared to the other two scenarios, is clear in Figure 11. It is evident that optimizing the locations of the monitoring wells regardless of existing observation wells in the study region leads to efficient collection of accurate groundwater levels. Our results also show that taking into account the existing observation wells provides useful information for optimizing the entire monitoring network.

Figure 11

Pareto fronts for the three scenarios obtained with the first run.

Figure 11

Pareto fronts for the three scenarios obtained with the first run.

Close modal

Under the second scenario, the main constraint is the fact that existing monitoring wells are used in the optimized network. This caused a rise of the RMSE in the network associated with Scenario 2 compared with the RMSE associated with Scenario 1. Choosing 10 wells, for instance, the RMSE values of the first and second scenarios were 1.180 and 1.489, respectively. Pareto fronts for the first and second scenarios are shown in Figure 11, where it is evident that the fully optimized monitoring network under Scenario 1 (where all monitoring locations are optimized) exhibits a lower RMSE than the Pareto fronts for Scenario 2.

It is seen in Figure 11 that the RMSE values of the first and third scenarios for 25 wells equaled 0.972 and 1.070, respectively. This vicinity of the two fronts implies that when the number of wells in the monitoring network exceeds 18, the number of wells is sufficiently high that all areas of the aquifer are covered. This is affirmed by the closeness of the first and third fronts in Figure 10 for well numbers equal to 28, 29, and 30.

Proper characterization of groundwater conditions relies on well-designed groundwater monitoring networks. This work introduced multi-objective optimization of groundwater monitoring network design with the evolutionary algorithm NSGA-II under three scenarios. The novel methodology was applied to the Eshtehard plain aquifer, Iran. The first part of this paper created time series of groundwater level values over the entire area of the plain with kriging, and stored the time series into a comprehensive data base. The second part of this paper optimized the network of observation wells employing the NSGA-II, and optimized groundwater levels were obtained with the IDW. The objectives of the optimization problem were the minimization of the number of monitoring wells and the minimization of the RMSE.

Three scenarios for groundwater network design were herein considered. Under the first scenario, all the monitoring wells' locations were optimized over the entire aquifer with the constraint being that not more than 30 wells would be deployed in the monitoring network. The second scenario involved redesigning an existing monitoring network and choosing the best subset of 18 existing wells that satisfies the multi-objective criteria. The third scenario involved expanding the monitoring network beyond the existing 18 monitoring wells. The key findings of this study are: the most efficient and accurate groundwater monitoring network design approach is that achieved under Scenario 1, when all the monitoring locations are optimized, regardless of existing monitoring wells.

The redesign results for Scenario 2 indicate that it is necessary to remove several wells among the existing wells of the network to achieve an optimized groundwater monitoring network according to the objectives of the design approach. Scenario 3 provides the solution of how to expand an existing groundwater monitoring network by choosing the optimal new wells.

This study relied on time series of groundwater levels to design groundwater monitoring networks. Previous studies on monitoring networks have not applied time series of groundwater levels because of complexities that arise in handling temporal variability within the spatial analysis. The optimization algorithm employed in this paper considered the entire area of the aquifer in search of the best monitoring sites. In brief, this study presented a groundwater monitoring network design method that could search all the aquifer area to find the best monitoring sites employing long-term groundwater-level data.

Aboutalebi
,
M.
,
Bozorg-Haddad
,
O.
&
Loáiciga
,
H. A.
2015
Optimal monthly reservoir operation rules for hydropower generation derived with SVR-NSGAII
.
Journal of Water Resources Planning and Management
141
(
11
),
04015029
.
Aboutalebi
,
M.
,
Bozorg-Haddad
,
O.
&
Loáiciga
,
H. A.
2016a
Application of the SVR-NSGAII to hydrograph routing in open channels
.
Journal of Irrigation and Drainage Engineering
142
(
3
),
04015061
.
Aboutalebi
,
M.
,
Bozorg-Haddad
,
O.
&
Loáiciga
,
H. A.
2016b
Simulation of methyl tertiary butyl ether concentrations in river-reservoir systems using support vector regression
.
Journal of Irrigation and Drainage Engineering
142
(
6
),
04016015
.
Aboutalebi
,
M.
,
Bozorg-Haddad
,
O.
&
Loáiciga
,
H. A.
2016c
Multi-objective design of water-quality monitoring networks in river-reservoir systems
.
Journal of Environmental Engineering
143
(
1
),
04016070
.
Akbari-Alashti
,
H.
,
Bozorg-Haddad
,
O.
,
Fallah-Mehdipour
,
E.
&
Mariño
,
M. A.
2014
Multi-reservoir real-time operation rules: a new genetic programming approach
.
Proceedings of the Institution of Civil Engineers: Water Management
167
(
10
),
561
576
.
Asefa
,
T.
,
Kemblowski
,
M. W.
,
Urroz
,
G.
,
McKee
,
M.
&
Khalil
,
A.
2004
Support vectors-based groundwater head observation networks design
.
Water Resources Research
40
(
11
).
DOI: 10.1029/2004WR003304
.
Bivand
,
R. S.
,
Pebesma
,
E.
&
Gómez-Rubio
,
V.
2008
Applied Spatial Data Analysis with R
.
Springer
,
New York
.
Bozorg-Haddad
,
O.
,
Soleimani
,
S.
&
Loáiciga
,
H. A.
2017
Modeling water-quality parameters using genetic algorithm-least squares support vector regression and genetic programming
.
Journal of Environmental Engineering
143
(
7
).
Cressie
,
N. A. C.
1991
Statistics for Spatial Data
.
John Wiley & Sons
,
Hoboken, NJ
.
Deb
,
K.
,
Pratap
,
A.
,
Agarwal
,
S.
&
Meyarivan
,
T.
2002
A fast and elitist multiobjective genetic algorithm: NSGA-II
.
IEEE Transactions on Evolutionary Computation
6
(
2
),
182
197
.
Dhar
,
A.
&
Patil
,
R. S.
2012
Multiobjective design of groundwater monitoring network under epistemic uncertainty
.
Water Resources Management
26
(
7
),
1809
1825
.
Dokou
,
Z.
&
Pinder
,
G.
2009
Optimal search strategy for the definition of a DNAPL source
.
Journal of Hydrology
376
(
3–4
),
542
556
.
Fallah-Mehdipour
,
E.
,
Bozorg-Haddad
,
O.
&
Mariño
,
M. A.
2013a
Extraction of optimal operation rules in aquifer-dam system: a genetic programming approach
.
Journal of Irrigation and Drainage Engineering
139
(
10
),
872
879
.
Fallah-Mehdipour
,
E.
,
Bozorg-Haddad
,
O.
,
Orouji
,
H.
&
Mariño
,
M. A.
2013b
Application of genetic programming in stage hydrograph routing of open channels
.
Water Resources Management
27
(
9
),
3261
3272
.
Fallah-Mehdipour
,
E.
,
Bozorg-Haddad
,
O.
&
Mariño
,
M. A.
2014
Genetic programming in groundwater modeling
.
Journal of Hydrologic Engineering
19
(
12
),
04014031
.
Loáiciga
,
H.
,
Charbeneau
,
R. J.
,
Everett
,
L. G.
,
Fogg
,
G. E.
,
Hobbs
,
B. F.
&
Rouhani
,
S.
1992
Review of ground-water quality monitoring network design
.
Journal of Hydraulic Engineering
118
(
1
),
11
37
.
Loáiciga
,
H. A.
,
Mariño
,
M. A.
&
Yeh
,
W.-G.
2010
Standard Guideline for the Geostatistical Estimation and Block-Averaging of Homogeneous and Isotropic Saturated Hydraulic Conductivit
.
ANSI/ASCE/EWRI Standard 54-10
,
ASCE Press
,
Reston, VA
.
McKinney
,
D.
&
Loucks
,
D. P.
1992
Network design for predicting groundwater contamination
.
Water Resources Research
28
(
1
),
133
147
.
Mogheir
,
Y.
,
de Lima
,
J. L. M. P.
&
Singh
,
V. P.
2003
Assessment of spatial structure of groundwater quality variables based on the entropy theory
.
Hydrology and Earth System Sciences
7
(
5
),
707
721
.
Orouji
,
H.
,
Bozorg-Haddad
,
O.
,
Fallah-Mehdipour
,
E.
&
Mariño
,
M. A.
2013
Modeling of water quality parameters using data-driven models
.
Journal of Environmental Engineering
139
(
7
),
947
957
.
Orouji
,
H.
,
Bozorg-Haddad
,
O.
,
Fallah-Mehdipour
,
E.
&
Mariño
,
M. A.
2014
Flood routing in branched river by genetic programming
.
Proceedings of the Institution of Civil Engineers: Water Management
167
(
2
),
115
123
.
Soleimani
,
S.
,
Bozorg-Haddad
,
O.
,
Saadatpour
,
M.
&
Loáiciga
,
H. A.
2016a
Optimal selective withdrawal rules using a coupled data mining model and genetic algorithm
.
Journal of Water Resources Planning and Management
142
(
12
),
04016064
.
Soleimani
,
S.
,
Bozorg-Haddad
,
O.
&
Loáiciga
,
H. A.
2016b
Reservoir operation rules with uncertainties in reservoir inflow and agricultural demand derived with stochastic dynamic programming
.
Journal of Irrigation and Drainage Engineering
142
(
11
),
04016046
.
Srinivas
,
N.
&
Deb
,
K.
1994
Multi-objective optimization using nondominated sorting in genetic algorithms
.
Evolutionary Computation
2
(
3
),
221
248
.
Yang
,
F.
,
Cao
,
S.
,
Liu
,
X.
&
Yang
,
K.
2008
Design of groundwater level monitoring network with ordinary kriging
.
Journal of Hydrodynamics
20
(
3
),
339
346
.