Abstract
This is the first attempt to combine the Multiobjective Billiards-Inspired Optimization Algorithm (MOBOA) with groundwater modelling to determine pumping rates within a well-distributed range of Pareto options. In this study, in order to determine an optimum solution for groundwater drawdown, pumping rates were selected accompanied by three minimization objectives: minimizing shortage influenced by inability to supply, adjusted shortage index, and minimizing the degree of drawdown within predefined areas. To optimize hydraulic conductivity and specific yield parameters of a modular three-dimensional finite-difference (MODFLOW) groundwater model, the Harris Hawks optimization algorithm was used to minimize the sum of absolute deviation between observed and simulated water-table levels. MOBOA was then utilized to optimize pumping rate variables for an Iranian arid to semiarid groundwater environment using these parameters. As the study results, when the maximum and minimum aquifer drawdown was specified in the range of −40 to +40 cm/year, the Pareto parameter sets produced satisfactory results. Overall, the “Simulation-Optimization-Modelling” protocol was able to generate a series of optimal solutions that were shown on a Pareto front. The study concluded to an optimum approach that provides policy makers in the Iranian water stressed zones with safe groundwater management alternatives.
1. Introduction
Groundwater behaviour modelling is one of the critical mechanisms that hydrogeologists have been attempting to measure since long time ago in order to solve evolving groundwater issues [1]. Due to the dynamic and multiobjective nature of the groundwater system, simulation of the groundwater system is challenging, especially in arid to semiarid zones. In the issue of rising groundwater demand, new models are desperately needed to develop novel decision-making tools and improve aquifer system drawdown [2]. Groundwater simulations have traditionally been conducted with the use of simulation/optimization algorithms [3, 4]. These models have been used to address construction and process issues of groundwater hydraulic control, water supply, and remediation [5, 6]. Professional experience with groundwater mechanism calibration shows which the single-objective functions are often insufficient to accurately quantify all dimensions and features of a groundwater system. Since sustainability managing is inherently a multiobjective issue, no optimum solutions can be determined in the conventional context, and policy makers can articulate their favourites through a collection of nondominated solutions [7].
To consider the multiobjective existence of groundwater systems, one approach is to specify multiple optimization objective functions that quantify different characteristics of system action. They employ a multiobjective optimization tool to find a collection of nondominated solutions named as Pareto best approaches [8, 9]. The Pareto solutions reflect trade-offs among various incomputable and often competing goals, with the property that switching from one solution to another improves one while deteriorating one or more others [10]. In groundwater research, there are many literature studies relating to the use of either deterministic or stochastic optimization approaches. Linear programming, nonlinear programming, and dynamic programming are examples of deterministic optimization approaches. Several academics have followed these approaches [11, 12]. The standard part of optimization techniques, known as metaheuristic algorithms, comprises genetic algorithms (GA), particle swarm optimization (PSO), simulating annealing, and others. These approaches were applied to extremely nonlinear or multimodal issues with a variety of dynamic constraints in order to improve groundwater characteristics in various regions [13].
Many attempts have been made in the literature to address multiobjective optimization of groundwater management issues. Park and Aral [14] proposed a multiobjective optimization method for determining coastal well positions that enhances pumping rates while minimizing the distance between the sensitive stagnation point and the reference coastline position. Reed et al. [15] developed a multiobjective strategy to cost-effective long-term groundwater monitoring using an Elitist nondominated sorted GA. In Italy, Giustolisi and Simeone [5] devised a method for assessing the complicated relationship between precipitation and water level in shallow unconfined aquifers. Siegfried et al. [16] proposed a multiobjective algorithm for optimizing pumping facility positioning and operation over time.
Saafan et al. [17] optimized pumping rates in Egypt’s El-Farafra oasis by combining a multiobjective genetic algorithm optimization model with MODFLOW. They forecast maximum pumping rate and minimum operation costs, as well as potential improvements in both variables. The above algorithms have concentrated on using a single genetic algorithm technique to calculate groundwater properties, which cannot be appropriate for large-scale groundwater systems. While working with a large-scale and complex structure, multiple competing priorities can occur, resulting in a rapid increase in the number of decision variables based on the problem scale. In such cases, single-objective approaches can produce unsatisfactory results for decision makers, necessitating the search for multiple optimization solutions. Many optimization domains include multidimensional problems that can only be solved concurrently by constructing different optimization algorithms.
El-Ghandour and Elsaid [12] suggested a steady-state analytical solution in a homogeneous unconfined aquifer using the particle swarm optimization (PSO) method to optimize and solve the groundwater management problem. They tested the proposed model on a popular hypothetical example to maximize the total pumping rate from located well system at steady-state condition. The results showed the superiority of the proposed model to obtain the maximum pumping rate compared with other methods of previous work.
It is the first research to combine a new multiobjective optimization algorithm called as Multiobjective Billiards-Inspired Optimization Algorithm (MOBOA) developed by Kaveh et al. [18] with MODFLOW to optimize pumping rates within a well-distributed range of Pareto approaches. The aim of the study was to determine the optimum approaches to meet the highest demand for arid to semiarid groundwater managing. Three minimization objectives were employed to optimize the pumping rates: minimizing the amount of drawdown within predefined Pareto zones, decreasing the shortage caused by an inability to supply, and decreasing the Modified Shortage Index (MSI).
Sadeghi-Tabas et al. [19] tried to link the multialgorithm genetically adaptive search method (AMALGAM) with a groundwater model to define pumping rates within a well-distributed set of Pareto solutions. The pumping rates along with three minimization objectives, i.e., minimizing shortage affected by the failure to supply, modified shortage index, and minimization of extent of drawdown within prespecified regions, were chosen to define an optimal solution for groundwater drawdown and subsidence in an arid groundwater system, Iran. The result was that the “Modeling-Optimization-Simulation” procedure was capable to compute a set of optimal solutions displayed on a Pareto front. In addition, the proposed optimal solution provides sustainable groundwater management alternatives to decision makers in arid region.
In addition, combination of MODFLOW with an advanced swarm-intelligence-based algorithm as the Harris Hawks optimization algorithm (HHO) developed by Abdel-Basset et al. [20] to calculate aquifer hydrodynamic parameters using an automated search process was the first innovation issue of this study. Therefore, The HHO approach’s capability resulted in more reliable calculations in error variance. In addition, it assisted in characterizing groundwater processes at each zone. Therefore, the aquifer properties were optimized in a highly parameterized model by developing the HHO approach within the groundwater numerical model.
The second aspect of innovation developed in this research (i.e., as the second step) was to identify optimal Pareto solutions for groundwater drawdown using a novel multiobjective algorithms named as the Billiards-Inspired Multiobjective Optimization Algorithm (MOBOA) model.
The recent optimization algorithm application guaranteed the groundwater control over an arid to semiarid area in Northeast of Iran.
In Figure 1, the procedure used in this research is illustrated graphically.
2. Material and Methods
2.1. Geographical Location of the Research Field
The Gorgan Plain is located in the arid to semiarid zone of Golestan Province, Iran. The Gorgan Plain aquifer system is located in the latitude and longitude of 36° W 37′ to 37° W 27′ North and 53° W 51′ to 54° W 51′ East (Figure 2). It encompasses an area of around 4393 km2. The Gorgan Plain is categorized as an arid to semiarid zone by the DoMarton climatic classification, with mean annual rainfall and temperature of 254 millimeters and 19°C, respectively. According to the mean sea level, the maximum and minimum elevations are 150 and −26 meters, respectively. The plain’s slope is steep in the south and mild in the north [21] Figure 1 illustrates the Gorgan Plain aquifer’s geographical location.
The aquifer’s bedrock is mainly hard rock in the south and saline water in the middle and northern parts. The majority of the land in the north is young alluvial fans, and clay flat stretches from the south to the north. The study area is mountainous in the south, and deserts and fine Caspian Sea sediments cover the western and northern regions. Agriculture tends to have the largest demand for water from March to July, with a subsequent decrease in demand. Furthermore, throughout the summer, drinking demand increases moderately, while industrial demand remains nearly steady during the year.
2.2. Groundwater Modelling
The mathematical solution of the groundwater model solves the math form of the mass balance equation in one area and generates a generally continuous approach to the surrounds. Every parameter value in terms of surface, volume, or time is represented by a distinct portion of the mass balance equation. In summary, a groundwater simulation mathematical model consists of a set of numeric value for different indicators in the balance equation. In other words, the balance equation is constructed for a specific aquifer zone, but it is generalized to the entire field [21]. Moreover, the results of groundwater modelling can be influenced by numerous mistakes, such as those caused by the groundwater conceptual model, the estimated solution of groundwater factors, and unexplained interactions between numerous factors and features. In fact, establishing a groundwater computational model necessitates the collection of various parameters such as hydraulic conductivity (k), transmissivity, and storage coefficient (or specific yield (Sy)). Groundwater modelling can produce a range of errors and misfits due to the scarcity of hydrogeology data (particularly in groundwater stressed regions), the heterogeneity of parameters in space and time, inaccurate configuration of aquifer characteristics (location, type, number of layers, distribution, etc.), and the scaling impact of variables. Since the numerical model is based on a conceptual model, therefore, the groundwater conceptual models frequently simplify real-world hydrogeology settings improperly [22]. Physical stability in groundwater mass and energy flows, which is a continuity mechanism, is usually described to allow predictions at discrete time points. All of these variables and complexity contribute to simulation bias and error, posing a challenge to modelling the groundwater problems, especially in arid to semiarid regions as the most groundwater stressed areas.
2.3. Governing Flow Equations of Groundwater Modelling
In groundwater simulation, the governing equation takes the following general form [22]:
In which, kx, ky, and kz are hydraulic conductivity tensors and h, Ss, and R are pressure head, specific storage, and recharge or discharge (positive and negative) aquifer elements, respectively. The thickness of the saturated layer of an unconfined aquifer correlates with groundwater level. Dupuit [22] considered that the horizontal flow should be governed in the whole aquifer as well as proportionality of the hydraulic gradient to the slope of the free surface [23]. The equation is based on Dupuit’s assumption if the flow is two-dimensional and transient, and the continuity equation is as follows [24]:
In which Sy represents the specific yield. The MODFLOW software was introduced in 1984 and has since been updated and improved to simulate in both steady and unsteady circumstances. A distributed hydrogeology model was developed in this study through using numerical data needed for simulating a functional relationship between prediction and observation data in order to model steady and unsteady-state conditions in the Gorgan Plain aquifer under various drawdown and control conditions.
2.3.1. Mathematical Representation of Aquifer Models
Izady et al. [24] suggested the method applied to construct the groundwater mathematical model in this analysis. The groundwater model, according to their approach, consists of the six following phases:(a)Acquiring all possible data;(b)Verifying and configuring control observations;(c)Identifying the geometry of the aquifer;(d)Primary estimate of the hydrodynamic parameters of the aquifer;(e)Description of the recharge and discharge values of the aquifer; and(f)Integrating the effects of previous phases producing the final conceptual model.
A one-layered unconfined aquifer with thicknesses ranging from 5 m to 55 m was examined as a conceptual groundwater model for the Gorgan Plain groundwater modelling. Surface elevation data, well logs, well locations and measurements, a geology map, hydrography, and recharge data were used to build the groundwater model. To identify the plain boundary situation, topography and geology maps were first used. The single optimization approach was then used to approximate the spatially distributed hydrodynamic characteristics. The temporal discharge difference was calculated in this analysis by pumping 33999 wells located inside the aquifer boundary (Figure 3). In Figure 3, the pumping well drilled in and distributed on Gorgan Plain aquifer has been shown. It should be noticed that this pumping well discharges the aforementioned aquifer.
In this analysis, seven data coverages were integrated into the MODFLOW model to create a groundwater numerical model, consisting of aquifer boundary circumstances, piezometers, pumping wells, surface recharges, drainage data, hydraulic conductivity, and specific yield. MODFLOW’s boundary conditions are determined by using a constant head boundary, a head based flux (river, drain, and general-head boundaries), and a known flux (recharge, evapotranspiration, wells, and stream). Because of the differential hydraulic gradients across multiple areas, the Gorgan Plain aquifer functions as a transient aquifer. In addition, this aquifer is recharged/discharged at inflow/outflow limits by adjacent aquifers, resulting in conditional and temporary results. In this research, water front entry cells or grids were used to identify recharge and discharge sites at inflow and outflow boundaries. For all of those borders, a specific-head-boundary state was considered to remain constant at each cell number in the model. As illustrated in Figure 3, only 281 of the 1340 observational wells with validated data, as monthly observed water-table level data, were utilized for calibration (Figure 4).
The number of pumping wells was 33999 pumping wells that consists of agricultural, drinking water, and industrial pumping wells. Furthermore, the aquifer contains 2764 springs and 336 Qantas, which entered in the MODFLOW. MODFLOW was originally simulated in steady state to describe homogeneous zones before being used to construct the aquifer model. The aquifer was separated into 30 homogeneous zones to determine hydraulic conductivity and specific yield values. In order to develop and construct a groundwater model, numerous settings were given, including model network architecture, stress cycle and time step selection, determination of starting position, boundary condition, and shape and number of saturated zones. To solve partial differential equation, the aquifer was separated into 7208 cells or grids due to geological heterogeneity. The plain’s gridding was constructed with 60 rows and 140 columns such that each layer cell is squared to 250 m × 250 m (Figure 5).
Cells outside of the aquifer boundary were then given a zero code, indicating that they have no influence on the simulation analysis. While individual model cell features may be set, combining cells with the same attributes into homogeneous zones significantly enhanced and simplified the modelling procedure. Figure 4 shows the modelled structure of the Gorgan Plain aquifer used by the MODFLOW model.
In this research, the modelling of groundwater flow in the Gorgan Plain was carried out over a one-year cycle, from October 2018 to September 2019 as the calibration period and from October 2019 to September 2020 as the validation period. 12-month tension cycles with a ten-day time phase were defined and used in the modelling methodology (i.e., for each stress cycle, three time measures were taken into account). To complete the numerical model, absolute values of bed rock depths, topography, and beginning water level data were interpolated by using the Kriging method and assigned to the network’s cells.
2.4. Harris Hawks Optimization Algorithm as the Single-Objective Optimization Approach
Modeling the Harris Hawks hunting technique yielded the Harris Hawks optimization (HHO) algorithm [25], a novel nature-inspired method. The algorithm’s operative and exploratory stages involve looking for prey, making a sudden pounce, and attacking in a variety of ways. The hunt is conducted at random using two exploring approaches. In the first way, Harris Hawks perch on a position that takes into account the positions of other family members and prey, whereas in the second way, the hawks wait on selected tall trees. Both strategies may be represented with the same probability of q as follows [20]:where and imply hawks’ position vectors in the present and subsequent iteration, respectively. is mentioned to a random hawk picked from the population. is the rabbit position. , , , , and are the numbers that was arbitrarily produced. and are higher and lesser bounds to generate accidental positions of the hawks’ habitat. implies the average position of hawks in the population which can be capable of being computed as follows [20]:where implies the i-th position vector of each hawk at t-th iteration and i = 1, …, h. is mentioned as the number of Harris Hawks’ population. According to escape energy E or fleeing of the rabbit, the algorithm’s transition from exploring to operation is possible as follows [20]:where denotes the primary rabbit energy randomly generated in [−1, 1]. specifies the maximum number of iterations. When E ≥ 1, hawks look for more regions to explore the rabbit’s locations as the exploration phase; otherwise, the exploitation phase is activated. The description of the failure ( < 0.5) or success ( ≥ 0.5) of the rabbit escape is accomplished in the algorithm, when the probability is the same. Furthermore, a forceful (E < 0.5) or gentle (E ≥ 0.5) besiege will be accomplished by hawks according to the rabbit energy. The gentle besiege can be expressed as follows [20]:where implies the dissimilarity between hawk and the rabbit positions, and the draw of the accidental strength of rabbit jump is completed by employing an arbitrary amount of rand. The definition of the forceful besiege, on the other hand, is as follows [19]:
If ( < 0.5) and (E ≥ 0.5) and as the rabbit can effectively flee, gentle besiege with advance rapid dives is done. The best possible dive can be selected by the hawks.
2.5. Multiobjective Billiards-Inspired Optimization Algorithm
This section provides an explanation of the basic principle of this novel physics-based metaheuristic. The mechanism of the Billiards-Inspired Multiobjective Optimization Algorithm is then outlined in the following section. The physics’ natural laws and the Billiards game embedded in the clash of balls are the major genesis of the BOA. Vector algebra and conservation laws govern when the balls collide with other balls. The kinetic energies of balls are preserved during collisions in perfectly elastic collisions, in addition to the sum of all angular velocity. Therefore, the following ultimate velocities of balls after crash in perpendicular and parallel directions to the impact line will be calculated as follows [18]:where and are the velocities of the first and 2nd balls during the collision and and are their velocities after the colliding. In addition, the symbols ║ and ┴ symbolize parallel and perpendicular elements, respectively. The masses of the balls are represented by the parameters m1 and m2. The linking vector’s unit vector is often denoted by ║. It is worth noting that the velocities’ perpendicular components stay constant, so the forces are merely directed to the collision axis, resulting in the perpendicular components of the angular velocities being preserved for balls. Therefore, the aforementioned equations indicate that if the masses of the balls are identical, the balls just transfer the parallel component of their velocities. The balls’ velocity components during the crash are shown in Figure 6 [18]:
The final positions of bodies defined in this algorithm using the equations of kinematics in the event of continuous accelerating are as follows [18]:
2.6. Billiards-Inspired Multiobjective Optimization Algorithm
All of the solution approaches that comprise multiple decision variables are modelled as the multidimensional balls in this algorithm. These balls are the MOBOA’s searching factors, and every dimension indicates a variable. In summary, the procedure begins with a random distribution of balls, and after that, it transfers nondominated strategies into an outside archive. Any members of the archive are chosen as pockets in each iteration. Following that, the balls are split into two collections: regular and cue balls. Every cue ball strikes its targeted ball, causing it to travel into a pocket. Whenever cue balls collide with each other, collision rules and vector algebra take over determining the movement of collided balls as well as their final status. MOBOA’s phases are illustrated below, and its pseudocode is shown in Figure 7 [18].
2.6.1. Phase 1: Initialization
In the search space, the initial population of multidimensional balls is generated at random as follows [18]:
In which, is the ith ball initial value. The vectors of variables xmax and xmin represent the maximum and minimum permissible values. rand[0,1] is an arbitrary variable, ranged by [0, 1], as well as the population number is defined by 2N.
2.6.2. Phase 2: Evaluation of Population
The objective functions each ball [f (x) = f1(x), f2(x), …, fk(x)] T should be evaluated [18].
2.6.3. Phase 3: Use an External Archive
The archive’s first mission is to save the best nondominated approaches discovered till now. All existing nondominated solutions are added to the archive in each iteration, whereas dominated solutions are deleted. To keep the archive size limit, a gridding process is being contemplated, which divides the objective space into hypercubes. When the archive capacity exceeds the maximum size limit, certain approaches in the packed hypercube are removed and some new solutions are added. This procedure broadens the range of the estimated Pareto front. The repository’s 2nd job is to include any candidates that deserve to be pockets [18].
2.6.4. Phase 4: Identifying the Pockets
Any representatives of the archive are chosen as pockets by means of the roulette-wheel selection approach. Equation (12) defines the suggested chance for selecting a hypercube to randomly draw a pocket [18]:where is the greater-than-zero selection pressure and ni defines the amount of solutions in the ith hypercube. As a result, staying in a less populated region enhances the chances of getting a hypercube. It should be noted that the user determines the amount of pockets, and these pockets increase the algorithm’s exploiting potential.
2.6.5. Phase 5: Arranging the Balls
In each iteration, the solution candidates are ordered in ascending order using the maximum strategy. An approach having a lower maximum value of fitness is more suitable because it is found in sparse regions. For the jth solution, the maximum fitness value is determined as follows [18]:
Following that, the ordered population is split into two equivalent classes (N pairs): (1) the upper and (2) the lower half as the regular and cue balls, respectively. Every cue ball corresponds to its pair in the regular balls category. This clustering approach creates the ideal situation for cue balls to exploit good-positioning regular balls.
2.6.6. Phase 6: Assigning Pockets to Balls
The lengths and cut angles between the balls and the pockets are used to assign the pockets to each pair of balls. These factors are taken into account while calculating the regarded shot score, which is calculated as follows [18]:
The shot score between the ith pair of balls and the mth bag is defined by . and are the ith regular ball and ith cue ball positions, correspondingly. In addition, Pm denotes the location of the mth bag, and M denotes the number of pockets. And the cosine of the cut angle is determined by , and it is measured using the following equation:
In which the dot product is denoted by the symbol “.”; based on the measured scores, all pair of balls is assigned to the highest scoring pocket. Another choice is to use the roulette-wheel selection method to calculate the pockets for every pair of balls [18].
2.6.7. Phase 7: Population Updating
Finally, a collision occurs between cue balls and regular balls. Regular balls are then driven into pockets by cue balls. As indicated by the equations below [18], the new locations of regular balls following collisions are near their pockets:
In which and are the new and old values of the ith regular ball’s jth vector, respectively. The vector is the jth variable of the mth pocket belonging to the ith pair of balls. The error rate is denoted by EROR; in addition, the iteration number and maximum iteration number are denoted by and , respectively. rand[−EROR, EROR] is an arbitrary number between [−EROR, EROR] [18]. Determining the ultimate location of cue balls following a collision is dependent on certain prerequisites, such as calculating the velocity of regular and cue balls. The velocity of the regular balls following the impact is shown as follows [18]:
In which is the velocity of the ith normal body following the crash, and is the ith body’s displacement vector. It is worth noting that the symbol “ˆ” denotes a corresponding unit vector. The acceleration rate is defined by the parameter , which is set to one. Cue ball velocities are calculated as follows [18]:
In which and are the ith cue ball’s velocities during the crash, correspondingly. In conclusion, the changed positions of the cue balls are calculated as follows [18]:
In which the factor determines a reducing variable that controls the movement of cue balls and has an initial value in the range [0, 1]. Figure 8 depicts the ball updating process.
2.6.8. Phase 8: Breaking Free from Local Optima
An escape limit, such as EL inside (0, 1), is known to avoid trapping in local optima. EL is compared to rand for each updated ball that is an arbitrary number regularly spread within (0, 1). If rand EL, the updated ball’s random dimension is regenerated as [18]
2.6.9. Phase 9: Evaluate the Boundary Restrictions
The ultimate location of balls may be put outside of the permissible ranges during the method of updating the location of balls. The dimensions of the balls that have been violated must be regenerated in this situation [18].
2.6.10. Phase 10: Evaluating the Termination Criteria
After a certain number of iterations, the search process will be terminated. If the condition is not met, the procedure is restarted at phase 2 [18]. This study used 5,000 function evaluations of 50 population sizes for optimization reasons. Then, to provide a reliable and computationally effective method to multiobjective optimization problems, 12 monthly discharge parameters were employed. After about 3,000 simulation cycles, the model met the convergence criterion. Pumping rates, as well as three minimization objectives, namely, minimizing shortages caused by a failure to supply, MSI, and minimizing the drawdown amount within predefined areas, were afterward selected to identify the best result for groundwater drawdown.
2.7. Developing the MODFLOW-HHO-MOBOA Structure
The MODFLOW-HHO-MOBOA framework is established. A simulator-optimizer model in the MATLAB software was used to evaluate and calculate aquifers’ properties. The simulator model (i.e., MODFLOW) and the optimizer algorithms (HHO-MOBOA) should be coupled in order to calibrate the aquifer’s hydrodynamic parameters and calculate pumping rates, and therefore, simulator model and the optimizer algorithms were linked by developing the interface subfunction in MATLAB that contributes to the simulator-optimizer model. In this research, Theo Olsthoorn’s [26] MODFLOW simulator model was used throughout the MATLAB platform (mfLab program). The MF2005NWT (MODFLOW-2005, Newton-Raphson formulation) software, which includes the Newton-Raphson solver to improve the outcome of unconfined groundwater-flow problems, as well as the Upstream Weighting Package (UPW), was also employed in this investigation for numerical solutions. UPW was used to determine aquifers’ features governing flow movement between cells in the MF 20005-NWT and MF-OWHM (MODFLOW-2005, One Water Hydrologic Flow Model) methods.
The MOBOA algorithmic was run by introducing an initial population and setting the algorithm parameters. This population is spread through a set of calculation nodes, which run the simulation model and calculate the objective function of the points obtained. The optimization process will be accomplished according to the flowchart proposed until the stopping criteria are encountered (Figure 1).
2.8. Evaluation of Model Results
Mean error (ME), mean absolute error (MAE), root mean square error (RMSE), and normalized root mean square error (NRMSE; nondimensional variants of the RMSE) were all utilized as criteria throughout the calibration procedure. ME, MAE, and RMSE/NRMSE all indicate error in the units (or squared units) of the constituent of interest, which aids in data interpretation [27, 28, 29]. The following equations (equations (21)–(24)) were used to estimate error:
In which , , and denote the observed head, simulated head, and mean value for observed head, respectively, n denotes the number of observational wells, and m denotes the number of monthly time stages (i ranges from 1 to 11 months, while t ranges from 1 to 12 months). The optimization model’s decision variables are dependent on the number of the aquifer discharge volumes (33999 wells) over 12 months. The three objective functions employed in this analysis are mentioned below. Minimizing the amount of shortage was considered as the first objective function [30]:
In which, denotes the aquifer discharge volume in the time span t (MCM), denotes the demand volume in the time period t (MCM), and n determines the cumulative number months. An MSI (equation (26)) was selected as the second objective function because it simply minimizes the total of shortages over time and avoids the spread of shortages over time. This index is critical for developing economic and social strategies [31]:
In which, TSt and TDt represent the amount of shortage and demand in period t, respectively. In addition, n indicates the total number months. It is worth noting that drinking demand is calculated by subtracting the shortfall volume of other needs (such as agricultural and industrial) from discharge values. Finally, the third objective function was to minimize the drawdown of water-table level [32]:
In which, and are the mean aquifer water levels (meters) at the start and end of the simulation time, correspondingly.
3. Discussion and Results
3.1. Optimizing Parameters Using the HHO Algorithm
Over 5,000 simulation runs, the HHO model was optimized, with the model achieving the convergence condition after around 3,000 iterations. For each zone, the specific yield and hydraulic conductivity of groundwater were then optimized. Figures 9 and 10 depict optimal hydraulic conductivity (k) and specific yield (Sy) values for geological units. The aquifer’s hydraulic conductivity and specific yield varies significantly, indicating a geologically nonhomogeneous groundwater system. The corresponding hydraulic conductivity field is greatest in blue-coloured zone and lowest in red-coloured one (Figure 11). It means that the central to the north western parts of the aquifer have the lowest hydraulic conductivity. Only two spots (i.e., blue spots) have the highest amount of hydraulic conductivity. In other words, the eastern parts of the Gorgan Plain aquifer have higher amount of hydraulic conductivity. As a matter of fact, this is due to the heterogeneity spread all over the aquifer. Although, due to the sever nonhomogeneity of the Gorgan Plain aquifer, the interpretation of hydraulic conductivity all over the study area is difficult, discussion about the specific yield can be more convenient. As it is obvious from Figure 9, the specific yield amount reduces to the northwest and northeast. Its highest amount located in the south with the blue colour as the blue spot. It means that the spreading pattern of specific yield would be concordant to the aquifer’s hydraulic conductivity except for the blue and green spots toward the east. Therefore, the before mentioned heterogeneities can be the main cause of this variation pattern.
High hydraulic conductivity regions can be permanently limited to the groundwater regime to modify aquifer properties. Considering the broad range of variations in both parameters across the entire aquifer, it seems that the aquifer’s hydraulic reaction to recharging and pumping operations is extremely volatile. To put it another way, there cannot be a steady-state pattern in the groundwater regime. The vertical range of the aquifer can be further disrupted by variations in hydraulic conductivity (hydraulic conductivity multiplied by saturated aquifer thickness). Hydraulic characteristics or other elements that have a significant positive or negative influence on groundwater modelling are of higher importance in the modelling [8]. The specific yield ranges between 0.02 and 0.2. In general, variations in specific yields can cause fluctuations in the aquifer’s hydraulic response to transient stress. This instability represents the groundwater regime’s extremely nonlinear characteristics.
After simulating each time, the mfLab software was slightly changed to display the volumetric specific yield. The specific yield values were separated by the research area to get equal water thicknesses. A variety of MATLAB commands were developed in this study to facilitate algorithmic runs and produce MODFLOW input files. All of the tests were performed on a PC using an Intel Core i7, 7700HQ CPU running at 2.8 GHz to 2.81 GHz and 12 GB of RAM. Depending on the algorithm setup, the overall running time was about 37 minutes that was much faster than using the internal optimization method supplied in MODFLOW software (i.e., Newton-Raphson optimization subfunction).
3.2. HHO Calibration and Validation Results
During the calibration process, aquifer water-table level was simulated across 30 zones, as seen in Figures 9 and 10. Using 12-month tension data from October 2019 to September 2020, the model was validated by predicting observed water-table level (each cycle includes three ten-day time scales). The contribution of different locations within the model domain to the estimated value of the particular parameter being measured is represented by the estimated values in distinct zones for a specific groundwater parameter. Table 1 shows and reports the error indices for simulated and observed head values. The error calculation indicates that the model is well calibrated. The calculated error is also low during the validation process, suggesting a good calibration. The MODFLOW simulation’s validity and the HHO’s efficiency in computing groundwater properties are also verified by the model calibration and validation results. Any ambiguity, such as measurement error, specified initial and boundary conditions, or the overall conceptual model, has a significant influence on modelling. Although Hamraz et al. [3] examined the parameter uncertainty associated in identifying the recharging ground water parameters; the emphasis of our study is on optimization approaches (by reducing error and misfit). According to the error indices, the HHO significantly decreased error and misfit by indicating the lowest error throughout the calibration and validation phases as well as reducing the CPU time. The HHO algorithm yielded RMSE and NRMSE in the ranges of 0.7–0.9 m and 0.02–0.025, respectively (Table 1).
The HHO approach discussed here will afford modellers with a simple technique that can be applied to complicated activities. As shown in Figures 10 and 12, the groundwater water-table levels were calculated at the start and end of the simulation phase for the whole aquifer.
The head values reveal a distinct pattern when comparing groundwater water-table levels. The effects of local pumping operations (drinking or agricultural wells) at the start of the simulation phase (which starts in October) and the subsequent restoration of flow can be used to explain this phenomenon (mostly by neighbouring aquifers). These aggressive pumping operations might smooth out the design of the simulated head curve at the aquifer’s centre. In other words, the values at the cell or grid’s centre are reflected in the simulated groundwater level, and piezometers are frequently found there. The aquifer’s groundwater level decreased from south to north and northwest, as seen in Figures 12 and 13, indicating that groundwater flow is largely in the same route. We can see from Figures 12 and 13 that if present withdrawals continue, groundwater levels in the west would almost certainly decline. In addition, the groundwater gradient decreases from south to north (and east to west in the eastern parts). This is due to the flow direction which is concordant to the groundwater gradient.
3.3. Pumping Rate Optimization
The optimized parameters collected by the HHO optimization algorithm were used to optimize predefined objective functions in MOBOA (5,000 iterations). Three policy analysis possibilities are examined to analyze the influence of different fundamental assumptions on the results. Scenario A: under this situation, the highest pumping rate is considered with reducing the sum of shortages while having the aquifer water-table drops maximally. Scenario B: in this situation, pumping rates are permitted to supply demand with no drop in aquifer water-table level and no volatility in the aquifer water table at the start and end of the year (October to September). Scenario C: pumping rates are expected to be the minimum values in order to minimize drawdown and raise the aquifer water level at the start and end of the period. Therefore, the MOBOA was optimized using 12-month discharge values as decision parameters according to the predetermined objective functions (Scenarios A, B, and C) as the alternative optimum approaches situated in separate areas of the Pareto solutions.
So, the MODFLOW framework is used to simulate and explain optimum solutions. Figure 13 depicts the Pareto-optimal solution using the three predefined objective functions values (i.e., drawdown, shortage, and MSI). As the distribution of shortages in different months is defined by MSI values, it (i.e., MSI) is considered to decline while shortage amounts reduce, and the drawdown amounts increase. Therefore, as drawdown increases, the shortage amount decreases, until all needs are supplied with no shortage when drawdown reaches 0.4 meter. Furthermore, 2D Pareto fronts (two objective functions) were evaluated to a better explanation and to analyze the best solutions (Figure 14). Drawdown and shortage are inversely related; because of the decreasing trend in the aquifer, discharge causes a drop in drawdown, as seen in Figure 14(a). The relation between shortage and drawdown is nearly linear, and it might be explained by the same amount for aquifer hydrodynamic characteristics at different flow depths. MSI has a straight association with shortage. As a result, MSI is projected to be inversely proportional to drawdown amounts (Figure 14(b)).
(a)
(b)
(c)
The total of optimum discharge amounts across the simulation time is shown as the decision variables for the three samples of optimum solution from various regions of the Pareto front. Therefore, the black, green, and blue dots in Figures 13 and 14, respectively, show the solutions A, B, and C.
According to Figures 13 and 14, if the first scenario (A) is considered, the total amount of shortages will be reduced and the aquifer will experience high values of drawdown. It will be accompanied by the total amount of MSI reduction as the highest pumping rate applied to the aquifer. It is well-known if Scenario C is considered, the aquifer water level rises upward due to declination in drawdown amount. This situation will not be obtained unless the aquifer pumping rates are expected to be at the minimum values.
Based on Scenario B, there will be no drop and instability in aquifer water level due to permitting the pumping rates to supply demands at the start and end of the year.
Figure 14 determines the solutions A and C that express zero, and 27.3 MCM shortages are decided to supply the first and third objective functions so as to minimize the sum of shortages and drawdown, respectively. The best solution would be one that all three objective functions contribute equally to the optimization procedure. The ideal solution is represented by the green point as solution B, when shortage and MSI are nearly zero. According to Figure 14, when drawdown is near 0.4 m, the aquifer can meet all demands. Moreover, in Scenario B, the pumping rate at existing wells is allowed to fluctuate over time (but not the water-table depth at the start and end of the period), resulting in a more arrangement that can respond to possible changes such as water needs or periodic recharges. Every value in the nondominated collections acquired for all options represents a distinct strategy that must be weighed against the others when making management choices. These options contrast from groundwater management viewpoints, which hold that no improvement can be made to one goal without decreasing satisfaction with the others. Figure 15 shows the amount of water-table level variations for the three chosen solutions across the simulation period. During February, these three best solutions are fairly near, but they are far apart in September.
February is an appropriate month since demand and drawdown are almost equal or have minimal variations. Although, throughout the spring and summer, once the aquifer water level drops steadily according to numerous pumping operations, the aquifer recovers significantly from October to February. The aquifer’s water-table level ranges from 32.6 to 32.8 meters in February but drops to 30.5 meters in September as the aquifer maximum discharge. During various options, the minimum and maximum limits of groundwater levels were 30.5 m and 32.9 m, respectively, according to the evaluation of all optimum solutions (Figure 15). As decision makers have three ideal alternatives for the Gorgan Plain aquifer based on the findings, they can choose any (drawdown) option that is within the ideal range of −40 to +40 cm/year. During maximum and minimum drawdown amounts over the aquifer, there are considerable changes in water-table levels. If one of the optimum management strategies, A, B, or C, is selected, it should be assured that the area’s maximum needs are satisfied. However, the challenge of deciding the best option stays open and may vary over time, based on management policies. It is worth noting that optimum solution B, as a transitional scenario, performed well in terms of decreasing the objective functions under consideration in this research. It demonstrates that the created technique was able to converge on an intermediate solution for the case study. Furthermore, optimal solution B could meet the supplying demand fairly while also reducing decline in the aquifer system.
The Pareto front’s optimal solution performances suggest that the assistances of the 281 wells might meet the requirements and quality standards. The best contributions produced the highest (+40 centimetres) and lowest (−40 centimetres) drawdown in the 281 wells at the end of the optimization phase. Figure 16 determines that the agricultural demand for water is the primary cause of water level fluctuations in the Gorgan Plain aquifer. In addition, the summer agricultural demand for water has been highlighted as the main discharge water.
It should be mentioned that the results of this research are in line with the Sadeghi-Tabas et al. [19] research. In other words, the results obtained in this study are similar to their study results.
4. Conclusion
In this research, the HHO method was used to optimize the estimation of MODFLOW hydrodynamic parameters, and the optimized parameters were then utilized in the MOBOA multiobjective algorithm to construct Pareto-optimum approaches and groundwater-managing scenarios. The suggested multiobjective problem design included reducing the shortage caused by a supply failure, MSI, and reducing the amount of drawdown within predefined zones within three alternatives as Scenarios A, B, and C. Under Scenario A, the total amount of shortages will be reduced and the aquifer will experience high values of drawdown. In addition, it will be accompanied by the total amount of MSI reduction as the highest pumping rate applied to the aquifer. However, if Scenario C is considered, the aquifer water level rises upward due to declination in drawdown amount. This situation will not be obtained unless the aquifer pumping rates are expected to be at the minimum values. Moreover, in Scenario B, there will be no drop and instability in aquifer water level due to permitting the pumping rates to supply demands at the start and end of the year. The results determined that considering the Scenarios A and C expresses zero and 27.3 MCM shortages, in order to supply the first and third objective functions so as to minimize the sum of shortages and drawdown, respectively. And the best solution would be one that all three objective functions contribute equally to the optimization procedure. However, the ideal solution is represented by Scenario B, when shortage and MSI are nearly zero. It is worth noting that, when drawdown is near 0.4 m, the aquifer can meet all demands. Moreover, in Scenario B, the pumping rate at existing wells is allowed to fluctuate over time (but not the water-table depth at the start and end of the period), resulting in a more arrangement that can respond to possible changes such as water needs or periodic recharges. In addition, the groundwater level variations for the three scenarios across the simulation period indicated that during February, these three best solutions are fairly near, but they are far apart in September. The aquifer’s water-table level also ranges from 32.6 to 32.8 meters in February but drops to 30.5 meters in September as the aquifer maximum discharge. During various options, the minimum and maximum limits of groundwater levels were 30.5 m and 32.9 m, respectively, according to the evaluation of all optimum solutions. When two separate algorithms were combined with the MODFLOW model, it was discovered that the created method could produce a set of optimum approaches for the groundwater system represented on a Pareto front. It should be noted that using the HHO algorithm as a new optimization method instead of the traditional optimization method supplied in the MODFLOW model makes the MODFLOW software more applicable and much faster in groundwater modelling by shortening the CPU time of the computers up to 37 minutes running time. Therefore, this research calculated optimal solutions for the whole aquifer of the Gorgan Plain, evaluating those approaches for each groundwater zone might give important data for the groundwater-managing policies. Although both the HHO and MOBOA algorithms are effective in quantifying groundwater features, the computation time for a broader aquifer system may grow considerably as the number of cells and stress periods rise. Therefore, the idea of fuzzy set theory [33] or using multiprocessor computers [34] can also be utilized to evaluate each of the efficient approaches to optimize objective functions depending on objective conditions. The HHO-MOBOA method described in this research can be regarded as an optimum solution producer, and its connection to the MODFLOW framework allows for the simulation of optimum groundwater scenarios that can be readily combined with the other conceptual models such as hydrology and water quality models. The integration of the optimization capabilities of HHO and MOBOA to the MODFLOW environment created a more applicable looping “Simulation-Optimization-Modelling” method valuable for decision makers dealing with groundwater management issues in the arid to semiarid regions.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request, e-mail address: [email protected].
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful to Dr. Alireza Donyaii for his kind recommendations. He is one of the best scientists in water resources managing in the Golestan Regional Water Company, Gorgan, Iran.