## ABSTRACT

This study aimed to determine the optimal conjunctive utilization of groundwater and surface water resources in the Halil River basin, one of the most significant study regions in Kerman Province (Iran). Multi-verse optimizer (MVO) and the ANFIS (adaptive neuro-fuzzy inference systems) simulation model, known as the MVO–ANFIS simulation–optimization model, were used for this purpose. Moreover, the optimal exploitation policy for the studied basin was presented. The ANFIS model yielded a coefficient of determination greater than 0.99 in Baft, Rabor, and Jiroft. This model had a high capability to simulate groundwater levels in these three regions. Therefore, the ANFIS model was adopted as the simulation model to predict the aquifer water table in these regions. Regarding the conjunctive utilization of groundwater and surface water resources, the exploitation policy resulting from the MVO–ANFIS simulation–optimization model had a desirable performance by supplying 91.70, 87.75, and 97.58% of the total demands of Baft, Rabor, and Jiroft, respectively. Moreover, results of water system performance indicators, including reliability (82.96, 72.65, 95.07), resiliency (70, 53.47, 80), vulnerability (29.54, 25.64, 17.02), and sustainability (74.24%, 66.10%, 85.78%) in the mentioned regions, respectively, showed the appropriate performance of the proposed model for the simulation–optimization problem.

## HIGHLIGHTS

An ANFIS model was used to simulate the underground water level.

The MVO–ANFIS model was developed for the purpose of conjunctive exploitation of surface and underground water resources.

The developed model had favorable results.

## INTRODUCTION

The lack of available water resources in Iran, a country with an arid and semi-arid environment, along with the growing water demand resulting from the expansion of agricultural, urban, and industrial activity, is consistently seen as one of the water sector's greatest concerns.

The water cycle refers to the continual flow of water between the earth and the atmosphere. Groundwater and surface water are part of the water cycle. Surface water involves any fresh water that flows into wetlands, rivers, and lakes. Groundwater is found in aquifers located underground. Most of the groundwater comes from snowmelt and rainfall and enters the bedrock through the surrounding soil. Groundwater and surface water resources are the main sources of water supply for human needs (Allan *et al.* 2020).

Comprehensive management of water resources includes the conjunctive utilization of groundwater and surface water resources. Preventing additional investment in surface water resources such as constructing dams or designing water transmission systems with an over-optimal capacity as well as excessive pressure on groundwater resources are among the benefits of conjunctive use. Conjunctive management of water resources refers to the conjunctive utilization of groundwater and surface water to optimize withdrawal from water resources (Zhang 2015).

Due to the complexity of the researched problem and objective functions, as well as the availability of multiple local and global optimal points, conventional approaches are unable to distinguish between them and locate the global optimal points in many engineering problems. Recent advancements in soft computing techniques such as genetic algorithm (GA) (Holland 1992), tabu search (TS) (Glover 1989), ant colony optimization (ACO) (Dorigo *et al.* 2006), particle swarm optimization (PSO) (Kennedy & Eberhart 1995), and multi-verse optimizer (MVO) (Mirjalili *et al.* 2016), which frequently search based on the initial population, have made it possible to find approximate solutions that are extremely close to the global optimum for extremely complex problems. Considering the capabilities of these methods and the necessity of conjunctive utilization of groundwater and surface water resources in regions with arid climates and shortages, using evolutionary methods will be useful.

Mahjoub *et al.* (2011) provided an integrated model of groundwater and surface water in the Maragheh Plain and used GMS (Groundwater Modeling System) to simulate 4 years of groundwater changes under drought and normal conditions. After modeling the groundwater of the Maragheh Plain, the pertinent data were incorporated into the integrated model. The groundwater and surface water supplies had no dynamic or hydraulic link. Milan *et al.* (2018) applied two fuzzy optimization techniques, including a fuzzy inference system and a linear fuzzy optimization model, to the conjunctive utilization of groundwater and surface water resources in the Astaneh-Kouchesfahan Plain. In the studied period, the mean shortfall values coming from the fuzzy inference system and the linear fuzzy optimization model were 22 and 14.6%, respectively, demonstrating the linear fuzzy optimization model's superiority. Under drought conditions, Seo *et al.* (2018) utilized a fully linked hydrological model for conjunctive management of groundwater and surface water resources. The results showed the proposed management model took a step towards the sustainable withdrawal of groundwater during drought periods.

Sepahvand *et al.* (2019) used a simulation–optimization model for conjunctive management of groundwater and surface water resources of the Gavkhuni basin considering wet, normal and dry planning periods. The results indicated the net profit of the agricultural sector increased in all three wet, normal and dry planning periods. Zeinali *et al.* (2020) investigated the conjunctive optimal use of groundwater and surface water resource systems by non-dominated sorting genetic algorithm II (NSGA-II). The results showed sustainable management of groundwater and surface water resources, especially in low-flow regions, taking into account various constraints.

Chakraei *et al.* (2021) investigated groundwater and surface water allocation using the NSGA-II multi-objective optimization model and multi-objective optimal planning approach for 10 years under different climatic conditions. Finally, the most appropriate solution was obtained from the interaction curve of objective functions of different stakeholders using the recursive bargaining model. Ashu & Lee (2021) used an artificial neural network (ANN) model integrated with Jaya algorithm (JA) optimization to determine the yearly conjunctive supply of agricultural water. The results showed this simulation and optimization model could minimize the water deficit by approximately 80%. Mani *et al.* (2016) employed the mixed integer linear fractional programming (MILFP) method for optimization use of groundwater and surface water. Results showed that the MILFP provides an efficient method to optimize conflicting objectives. Shukri *et al.* (2021) compared two algorithms, MVO and PSO3, using meta-heuristic algorithms to calculate task scheduling in computing environments. The results showed the MVO algorithm was significantly better in terms of achieving the minimum time to perform functional tasks. Jain *et al.* (2019) used the MVO algorithm for the automatic generation of computer programs. This study showed that the MVO algorithm can be used in the automatic generation of computer programs in any target computer language.

Recently, different novel hybrid machine learning algorithms were used in studies such as least square support vector machine, the improved version of multi-verse optimizer algorithm (LS-SVM-IMVO) (Adnan *et al.* 2022), Long short-term memory network- weIghted meaN oF vectOrs (LSTM-INFO) (Adnan *et al.* 2023a, 2023b), relevance vector machine tuned with improved Manta-Ray foraging optimization (RVM-IMRFO) (Adnan *et al.* 2023a, 2023b), adaptive neuro-fuzzy inference system, water cycle optimization algorithm and moth flame optimization algorithm (ANFIS-WCAMFO) (Adnan *et al.* 2021), extreme learning machines–jellyfish search optimizer (ELM–JFO) (Adnan *et al.* 2023a, 2023b), support vector machine, fruit fly optimization algorithm and particle swarm optimization (SVM–FFAPSO) (Adnan *et al.* 2022). In light of the growing need for maximizing the exploitation of water systems and the implementation of fresh and enhanced MVO algorithms in other fields of study, this study assessed the efficiency of MVO algorithms in optimizing the allocation of water resources. Moreover, ANFIS is a simulation model for underground water level and because it has been used before, we did not use it in our research, we used the ANFIS-MVO hybrid model and checked its results compared with ANFIS-GA and ANFIS-PSO methods because this method has not been used in water engineering. In the programming portion of the MATLAB software, an MVO-based model was constructed for this purpose. Conjunctive utilization is mostly used for places that do not have enough surface water and have to use ground water in a combined manner in arid and semi-arid areas and considering that the area we studied is a relatively semi-arid area, and all our needs cannot be completely met with surface and ground water sources, we have to go towards conjunctive utilization of surface and underground sources.

Using reliability, resiliency, and vulnerability as reservoir performance metrics, the created model's outcomes were evaluated.

## MATERIALS AND METHODS

### Study area

^{2}consists of 21 study areas and four large cities. In terms of the hydrological division of Iran, the Jazmurian basin is a part of the central desert basin of Iran and is divided into Lut and Daranjir-Saghand desert basins from the north, Hamun-Mashkel basin from the east, Abarqu-Sirjan desert basin from the west, river basin between Bandar Abbas and Sedij, Baluchistan rivers basin and Kol-Mehran basin. This is geographically located from 56° 15′ to 61° 23′ east longitude and 26° 28′ to 29° 30′ north latitude. In terms of country division, it includes parts of Kerman Sistan and Baluchistan provinces. Baft, Jiroft, Bampur and Iranshahr are located in this basin. Parts of sub-roads from Bam to Minab and Bam to Hajiabad as well as Khash to Chabahar road are located in the west and east of this basin, respectively. The southern part is limited. The important dams in the western part of this basin include the Baft earth dam (Asiab Jafteh) for the purposes of agriculture, drinking water and industry, the Jiroft double-arch concrete dam for the purposes of drinking water, industry, agriculture, electricity generation and flood control and Safaroud Dam (under construction) for the purposes of agriculture, drinking water and industry. Figure 1 depicts the position of the aforementioned dams in the western portion of the Jazmurian basin. (Halil River basin).

### Simulation model

#### ANFIS

Rule 1: If

Rule 2: if

The previously stated regulations can also be restated using the concept of membership degrees as follows:

Rule 1: If

Rule 2: If

• and are membership functions (MFs)

• , , are the parameters of the consequent part of fuzzy rules

The adaptable nodes in the ANFIS architecture are found in layer 1 (MFs layer) and layer 4 (consequent layer), while the fixed nodes are located in layer 2 (product layer) and layer 3 (normalization layer). Let us explore the five-layer structure of ANFIS in more detail:

*c*, ) serve as the premise parameters and are adjusted throughout the learning process.

Using grid partitioning, the automatically generated rules amount to mn, with ‘m’ representing the number of MFs for each input and ‘*n*’ being the total count of inputs.

- Layer 3: Every node, indicated by
*N*, normalizes a rule's firing strength by determining the proportion of the*i*th rule's firing strength relative to the total firing strength of all rules.where , and*w*correspond to the firing strength of the*i*th rule, the firing strength of the first rule, the firing strength of the second rule, and the normalized firing strength of the*i*th rule, respectively.

In the context given, *w* represents the normalized firing strength of the *i*th rule and represents a first-order polynomial of the consequent part of the *i*th rule. The parameters are determined through the training process of ANFIS

#### Developing ANFIS simulation model

Given that the equations regulating groundwater flow are based on the law of conservation of mass, it is necessary to establish the inputs and outputs driving water table changes. Effective inputs and outputs during monthly time steps in the study area are as follows:

Amount of groundwater withdrawal in the current time step (

*G*)._{t}Amount of allocation to the needs of each region in the current time step (De

)._{t}Precipitation rate in the current time step (Pr

)._{t}Mean monthly evaporation in the current time step (Ev

)._{t}Mean river flow (hydrometric station data) in the current time step (

*Q*)._{t}To include the initial conditions, the mean water level of piezometers of each region in the previous time step is considered as the input of the ANFIS simulation model (G-Lev

_{i}_{, t−1}).To add aquifer boundary conditions, the mean water level measured by piezometers in neighboring locations during the preceding time step is used as an input to the ANFIS model. In this way, the effect of changes in the groundwater level in each region on the conditions of the adjacent regions is taken into account (G-Lev

_{n}_{, t−1}).Finally, the ANFIS model output is the mean groundwater level of the region in the current time step (G-Lev

)._{t}

Therefore, there are seven inputs and one output for each region in the ANFIS model.

#### MVO

*et al.*2016) that is based on the multi-verse theory, which suggests that the universe was generated by several big bangs. This theory explains the existence of several parallel universes. The MVO algorithm is based on three primary components: the white hole, the black hole, and the wormhole (Figure 3).

#### MVO algorithm

The following principles are applied to MVO universes during optimization:

- 1.
The greater the inflation rate, the greater the likelihood of white holes.

- 2.
The higher the rate of inflation, the lesser the likelihood of black holes.

- 3.
In universes with a higher inflation rate, items are more likely to travel through white holes.

- 4.
Universes with lower inflation rates typically receive a greater number of things through black holes.

- 5.
Regardless of the inflation rate, objects in all universes can flow through wormholes to the universe with the best conditions.

*d*is the number of parameters (variables).*x*is the_{j}*j*th parameter in the*i*th universe.*u*is the_{i}*i*th universe.*NI*(*U*) is the normal inflation rate in the_{i}*i*th universe.*r*1 is a random number between [0, 1].*xk*,*j*is the*j*th parameter in the*k*th universe chosen by the roulette wheel selection method.

*x*is the_{j}*j*th parameter in the best universe formed to date.TDR and WEP are coefficients.

*lb*is the lower bound of the_{j}*j*th variable.*ub*is the upper bound of the_{j}*j*th variable.*r*2,*r*3, and*r*4 are random numbers between 0 and 1.

min and max represent the minimum (0.2 in this paper) and maximum (1 in this paper), respectively.

*l*is the current iteration.*p*(equal to 6 in this study) defines the accuracy of iterative exploitation.

*p*value, the more rapidly and precisely the local search and exploitation is accomplished.

*n*log

*n*) and O(

*n*2), respectively. Depending on the implementation, the roulette wheel selection is executed for each variable in each universe across iterations with an O(

*n*) or O(log) complexity. Therefore, the following describes the overall computational complexity:

*n*is the number of universes.*l*represents the maximum number of iterations.*d*represents the quantity of items.

The following factors are evaluated to determine whether the suggested algorithm has the theoretical capacity to solve optimization problems.

White holes are more likely to form in universes with a high rate of inflation. Thus, they may transport things to other worlds and assist them in increasing their rates of inflation.

In universes with modest rates of inflation, black holes are more likely to form. Thus, they have a greater chance of receiving things from other universes. This issue raises the likelihood of inflation rates rising in universes with low inflation.

Objects are typically transported through black/white hole tunnels from universes with high inflation rates to universes with low inflation rates. Thus, the total/mean inflation rate improves over iterations in all universes.

Wormholes form at random in all universes, independent of the rate of inflation. Thus, the variety of worlds might be preserved across iterations.

Black/white hole tunnels necessitate an abrupt transition in universes, which could ensure exploration of the search space.

Sudden changes overcome the problem of stagnant local optimums.

Wormholes randomly re-span some variables of universes around the optimal solution achieved by iterations, hence ensuring exploitation in the optimal search space.

Adaptive WEP values improve the likelihood that wormholes exist in the cosmos in a continuous manner. Consequently, the optimization approach emphasizes the exploitation process. Adaptive TDR values decrease the travelling distance of variables surrounding the best universes, hence enhancing the local search precision across iterations.

The convergence of the suggested method is ensured by stressing exploitation/local search in proportion to the number of iterations.

### General structure of integrated groundwater and surface water optimization model

After building the simulation model, the optimization model is created in order to achieve the customized optimal values for groundwater and surface water by combining the two models. The optimization model requires following the constraints applied to the state and decision variables. Thus, groundwater level fluctuations are considered as the state variable in the present study. The optimization algorithm produces decision variables and, then, transfers them to the simulation model obtained using the ANFIS model. The simulation model describes the system behavior in response to the made decisions and calculates the state variable value. The constraints are satisfied using the penalty method. If the generated solutions cannot satisfy the constraints and do not obtain the minimum objective function value, this process is repeated, so that the optimization algorithm regenerates new decision variables and sends them to the simulation model. This cycle continues until convergence to the global or near-global optimum.

### Optimization model objective function

,

*F*is the objective function.*i*is the time step counter (number of months in the planning period).*k*is the region counter.demand

is the net water demand in the_{k,i}*k*th region and*i*th month.supply

is the quantity of water allocated to the_{k,i}*k*th region in the*i*th month by the optimization model.*G*is the quantity of groundwater withdrawal in the_{k,i}*k*th region and*i*th month.*S*is the quantity of surface water withdrawal in the_{k,i}*k*th region and*i*th month.

### Constraint

*S*is the amount of monthly surface water withdrawal in the_{k,i}*k*th region and*i*th month.*S*_{max k,i}is the maximum amount of monthly surface water withdrawal in the*k*th region and*i*th month.*i*_{max,k}is the maximum amount of total surface water withdrawal in a planning period (223 months) in the*k*th region.

*L*is the groundwater depth, i.e., the distance between the well opening and aquifer level in the_{i,k}*i*th month and*k*th region.Δ

*L*is the change in the water level in the ith month compared to the (_{i,k}*i*− 1)th month in the*k*th region.Δ

*L*-total_{max,k}is the upper bound of the drop level or lower bound of the allowed water level improvement at the end of the planning period.

### Evaluating results using performance indicators

Planning policies and water resources management aim to reduce the impact of policies that negatively affect water resources systems both in the current situation and future and develop policies that have positive socioeconomic, environmental, political and legal impacts on systems. Therefore, there should be parameters or indicators for measuring system performance and evaluating and comparing water resources system conditions under different management policies and scenarios (Sandoval-Solis *et al.* 2011). Using optimization and simulation models to exploit a reservoir system culminates with the most crucial stage of evaluating exploitation policies. In this research, reliability, vulnerability, resiliency and sustainability indicators are used to evaluate the studied algorithms in the reservoir system.

### Reliability

NDe

is the total number of failures during the exploitation period._{f}De

is the demand value at time_{t}*t*, Reis the output value at time_{t}*t*and Rel is the system reliability during the exploitation period.

The bigger the value of this parameter, the stronger the system time reliability would be (Hashimoto *et al.* 1982).

### Vulnerability

*et al.*1982).

Val represents the system's vulnerability.

De

represents the demand value at time_{t}*t*.Re

represents the output value at time_{t}*t*.*T*represents the total number of exploitation periods.

### Resiliency

*et al.*1982). This parameter's value could be computed with Equation (22).

is the number of periods, in which the condition is provided in parentheses.

Def

is the deficiency at time_{t}*t*.

### Sustainability

*et al.*2011):

### Optimization model performance

After determining the initial parameters of the MVO–ANFIS simulation–optimization model, the integrated model is implemented. In the optimization algorithm, values are assigned to decision variables (withdrawal rate from groundwater and surface water resources) in all time-steps for each universe in the population. These values along with other inputs of the simulation model (which are done using the ANFIS model) are entered into the ANFIS model and the groundwater level, which is the output and state variable of the optimization model, is estimated. Since the optimization model should comply with the constraints applied to the state and decision variables, the constraint satisfaction is controlled using the penalty method. Finally, a value is given to the function of each universe whether the constraints are met or not. The universes are scored based on the obtained objective function values.

Due to the demand for water, this demand is supplied using the presented management model, which reduces water waste in some months and compensates for water scarcity in other months.

## RESULTS AND DISCUSSION

In this section, the results of executing the optimal exploitation model for groundwater and surface water resources in the Halil River basin (Baft, Rabor, and Jiroft) with the goal function of minimizing total demand insufficiency over a 223-month period are reported (2000–2019).

Table 1 displays the optimal objective function value, supply percentage, and demand deficiency in Baft, Rabor, and Jiroft after 1,000 iterations of the constructed MVO–ANFIS, PSO–ANFIS, and GA–ANFIS models on the objective function during the long period.

Method . | Objective function . | Baft . | Rabor . | Jiroft . | |||
---|---|---|---|---|---|---|---|

Deficiency (MCM) . | Supply percentage . | Deficiency (MCM) . | Supply percentage . | deficiency (MCM) . | Supply percentage . | ||

MVO-ANFIS | 794.14 | 19.58 | 91.70 | 69.42 | 87.75 | 28.50 | 97.58 |

PSO–ANFIS | 950.89 | 19.89 | 91.57 | 70.69 | 84.65 | 29.46 | 96.52 |

GA–ANFIS | 1,486.89 | 21.70 | 90.80 | 73.47 | 80.56 | 43.29 | 95.50 |

Method . | Objective function . | Baft . | Rabor . | Jiroft . | |||
---|---|---|---|---|---|---|---|

Deficiency (MCM) . | Supply percentage . | Deficiency (MCM) . | Supply percentage . | deficiency (MCM) . | Supply percentage . | ||

MVO-ANFIS | 794.14 | 19.58 | 91.70 | 69.42 | 87.75 | 28.50 | 97.58 |

PSO–ANFIS | 950.89 | 19.89 | 91.57 | 70.69 | 84.65 | 29.46 | 96.52 |

GA–ANFIS | 1,486.89 | 21.70 | 90.80 | 73.47 | 80.56 | 43.29 | 95.50 |

As presented in Table 1, the objective function values in the developed MVO–ANFIS, PSO–ANFIS, and GA–ANFIS models were 794.14, 950.89, and 1,486.89, respectively. This result indicates the high capability of the developed MVO–ANFIS model in solving high-dimensional problems compared to the other models. By jointly managing groundwater and surface water resources, the MVO–ANFIS technique met 91.70, 87.75, and 97.58% of Baft, Rabor, and Jiroft's total demands, respectively. For the next rank, the PSO–ANFIS, and GA–ANFIS models showed the demand-supply of (91.57, 84.65, 96.52%) and (90.80, 80.56, and 95.50%), respectively. Arya Azar *et al.* (2022) proposed FA-GMDH-LS-SVM and WOA-GMDH-LS-SVM models for simulation, optimization, and estimation of conjunctive use of surface water and groundwater resources. The water supply of these models was 87.1 and 84.4%. These results show that despite the fact that the models used in the mentioned study and the present study were used for different geographical areas, the MVO–ANFIS technique had better results. However, Zeinali *et al.* (2020) used the NSGA-II-WEAP-MODFLOW model with a water supply of 97.85%.

As indicated in Figure 7, the highest curve in the figure is represented by the black line with triangular markers, which is labeled as MVO–ANFIS. The lowest curve is depicted by the blue line with circular markers, representing the real-condition operation. This indicates that the MVO–ANFIS followed by PSO–ANFIS, and GA–ANFIS managed to maintain the groundwater level at a higher elevation compared to the actual condition. In simpler terms, the non-optimal condition caused the groundwater level to drop approximately 2 m more than the optimal condition in all the regions studied. This demonstrates that the optimal management scenarios were able to preserve a significant volume of groundwater. Consequently, it can be concluded that all the developed models successfully optimized the combined utilization of surface and groundwater resources in the Halil River basin, ensuring a reliable water supply. Furthermore, these scenarios effectively contributed to regulating the groundwater level and ensuring its sustainability in the study areas over a 19-year period.

Table 2 summarizes the evaluation criteria for the MVO–ANFIS, PSO–ANFIS, and GA–ANFIS simulation–optimization models for optimal conjunctive utilization of groundwater and surface water resources in the Halil River basin, including reliability, resiliency, vulnerability, and sustainability during the exploitation period.

Region . | Algorithm . | Reliability . | Resiliency . | Vulnerability . | Sustainability . |
---|---|---|---|---|---|

Baft | MVO-ANFIS | 82.96 | 70 | 29.54 | 74.24 |

PSO–ANFIS | 73.99 | 66.17 | 45.58 | 64.35 | |

GA–ANFIS | 73.09 | 65.71 | 66.01 | 54.65 | |

Rabor | MVO-ANFIS | 72.65 | 53.47 | 25.64 | 66.10 |

PSO–ANFIS | 69.05 | 57.77 | 45.42 | 60.16 | |

GA–ANFIS | 65.47 | 65.55 | 64.26 | 53.53 | |

Jiroft | MVO-ANFIS | 95.07 | 80 | 17.02 | 85.78 |

PSO–ANFIS | 85.65 | 61.36 | 24.79 | 73.39 | |

GA–ANFIS | 93.27 | 82.35 | 58.90 | 68.09 |

Region . | Algorithm . | Reliability . | Resiliency . | Vulnerability . | Sustainability . |
---|---|---|---|---|---|

Baft | MVO-ANFIS | 82.96 | 70 | 29.54 | 74.24 |

PSO–ANFIS | 73.99 | 66.17 | 45.58 | 64.35 | |

GA–ANFIS | 73.09 | 65.71 | 66.01 | 54.65 | |

Rabor | MVO-ANFIS | 72.65 | 53.47 | 25.64 | 66.10 |

PSO–ANFIS | 69.05 | 57.77 | 45.42 | 60.16 | |

GA–ANFIS | 65.47 | 65.55 | 64.26 | 53.53 | |

Jiroft | MVO-ANFIS | 95.07 | 80 | 17.02 | 85.78 |

PSO–ANFIS | 85.65 | 61.36 | 24.79 | 73.39 | |

GA–ANFIS | 93.27 | 82.35 | 58.90 | 68.09 |

As indicated in Table 2, the sustainability index which is a criterion for identifying the performance of water resource systems was obtained by MVO–ANFIS, PSO–ANFIS, and GA–ANFIS models as (74.24, 66.10 and 85.78%), (64.35, 60.16 and 73.39%), and (54.65, 53.53 and 68.09%) to supply demands of Bafat, Rabor and Jiraft, respectively. Therefore, it could be stated that the created MVO–ANFIS compared with the other models possessed a high level of capability and precision in solving the problem of optimal conjunctive utilization of groundwater and surface water resources to meet downstream needs. Akbarifard *et al.* (2024) employed the SOSMSA–ANN, MSA–ANN, and SOS–ANN models for the conjunctive operation of surface and groundwater resources in the same basin of this study. The results of their study showed that the sustainability index that was obtained by SOSMSA–ANN, MSA–ANN, and SOS–ANN models were (89.30, 66.69, and 57.59%), (81.58, 70.58, and 79.48%), and (97.06, 88.42, and 75.13%) to supply demands of Bafat, Rabor and Jiraft, respectively. These results show that the models used in the mentioned study obtained a better level of sustainability than the models used in the present study.

## CONCLUSION

The ANFIS simulation model utilized in this research could predict the aquifer water level in different regions with an acceptable error rate, so that the determination coefficients (*R*^{2}) of Baft, Rabor and Jiroft were obtained as 0.9965, 0.9952 and 0.9996, respectively. Using the ANFIS model could be useful when it is not possible to integrate simulation and optimization models due to its high accuracy. Moreover, it would be easier to integrate the optimization models with the ANFIS model using the method presented in this research, which could increase the chance of approaching the global optimum solutions in a conjunctive management problem with linear and non-linear objective functions. The results obtained from implementing the MVO–ANFIS model showed although demands were fully supplied in some months, the proposed model could increase the mean aquifer water table to an acceptable level and provide acceptable performance indicators while satisfying all the constraints. The sustainability index of the hybrid MVO–ANFIS model was obtained as 74.24, 66.10, and 85.78% for supplying demands of Baft, Rabor and Jiroft, respectively, showing the model's suitability for the researched situation. The model could reduce the sharp increase and decrease of the aquifer water level in different months during the operation period, which could stabilize the aquifer water level and prevent land subsidence.

MVO performed well in reaching global or near-global optimal solutions. Moreover, the solution of not getting trapped in the local optimum was included in the MVO structure. Using MVO in solving non-linear problems (like the one investigated in this study) and examining its results help us whether to use MVO or not in solving complex problems that could not be solved by the existing optimization software. According to the investigations and obtained results, it is recommended to use the MVO algorithm and integrate it with the ANFIS optimization model for solving simulation–optimization problems. It is also recommended that researchers check the noise and uncertainty in the data to estimate the data in future studies.

## AUTHORS’ CONTRIBUTION

S.H. and A.R. contributed to designing the study. The data collection and analysis were done by S.H. S.H., A.R., M.M., and M.K. participated in drafting the manuscript. All authors read and approved the final version of the manuscript.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.