Prediction of lake water-level fluctuations using adaptive neuro-fuzzy inference system hybridized with metaheuristic optimization algorithms

Lakes help increase the sustainability of the natural environment and decrease food chain risk, agriculture, ecosystem services, and leisure recreational activities locally and globally. Reliable simulation of monthly lake water levels is still an ongoing demand for multiple environmental and hydro-informatics engineering applications. The current research aims to utilize newly developed hybrid data-intelligence models based on the ensemble adaptive neuro-fuzzy inference system (ANFIS) coupled with metaheuristics algorithms for lake water-level simulation by considering the effect of seasonality on Titicaca Lake water-level fluctuations. The classical ANFIS model was trained using three metaheuristics nature-inspired optimization algorithms, including the genetic algorithm (ANFIS-GA), particle swarm optimizer (ANFIS-PSO), and whale optimization algorithm (ANFIS-WOA). For determining the best set of the input variables, an evolutionary approach based on several lag months has been utilized prior to the lake water-level simulation process using the hybrid models. The proposed hybrid models were investigated for accurately simulating the monthly water levels at Titicaca Lake. The ANFIS-WOA model exhibited the best prediction performance for lake water-level pattern measurement in this study. For the best scenario (the inputs were Xt-1,Xt-2,Xt-3,Xt-4,Xt-12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X}_{t-1},\; {X}_{t-2}, \;{X}_{t-3}, \;{X}_{t-4}, \; {X}_{t-12}$$\end{document}) the ANFIS-WOA model attained root mean square error (RMSE ≈\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx$$\end{document} 0.08 m), mean absolute error (MAE ≈\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx$$\end{document} 0.06 m), and coefficient of determination (R2≈\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx$$\end{document} 0.96). Also, the results showed that long-term seasonal memory for this lake is suitable input for lake water-level models so that the long-term dynamic memory of 1-year time series for lake water-level data is the best input for estimating the water level of Titicaca Lake.


Introduction
Lakes were considered as one of the most important landscapes around the world, it is well known from one country to another that the hydrological conditions significantly influence the global functions of lakes (Ye et al. 2020). An in-depth understanding of lake water-level (LWL) fluctuations has been involved in the long-term planning and management of water resources. The direct link between climatic variables (e.g., precipitation and temperature) and LWL fluctuations must be accurately and adequately established (John and John 2019) to develop better models of LWL. LWL fluctuations were considered to be related and largely influenced by climate change for a long period, and it was observed that climate change influences precipitation, temperature, and evaporation regimes, which ultimately pushes the LWL to extremities (Munyasya et al. 2022, Yin et al. 2022, Bhele et al. 2020. The qualities, availability, complexity, and biological events among fauna are greatly influenced by the increasing water-level fluctuations (Gownaris et al. 2018). Therefore, effective management of lakes should be accompanied by proper monitoring of LWL fluctuations and an in-depth understanding of the effect of these fluctuations on the aquatic biome (e.g., benthic macro invertebrates) (Yan et al. 2020). In addition, LWL fluctuations also exercise direct and indirect effects on the physical, chemical, biological, and biogeochemical characteristics of lakes (Liu 1 3 13 Page 2 of 17 et al. 2019). Thus, the LWL regulation is very beneficial to assess flood risk and mitigate its impact . It is also imperative to monitor the fluctuation of LWL at different time scales, with particular emphasis on improving the forecasting methods, which intend to help develop early warning systems.
The accurate monthly LWL simulation remains up-todate on ongoing research interest for several water engineering disciplines, including water resources optimization, determining the permitted release rate of pollutants, and preparing for potential flood hazards and drought periods (Rezaie-Balf et al. 2017). Prediction of LWL at different time intervals using historical time series data is one of the challenging problems in water resource planning . Changes in LWL are affected by several factors such as weather variables (e.g., precipitation, temperature etc.), direct and indirect runoff from adjacent basins, evaporation from the free surface of the lake; however, others stem from the interactions between the lake and the groundwater table. Lake-level surveys have always been of interest to hydrologists and environmental experts. The volume of lake water, which impacts the management of groundwater aquifers, ecological and environmental changes in the region, is affected by the LWL. Although it is necessary to pay attention to the above parameters to identify complex models, it is better to use the LWL recorded data for monitoring its changes . The Titicaca Lake (TL) water-level fluctuations significantly impact water resource planning and management, freshwater supply management, fishing, and other recreational activities (Kroll et al. 2012) in both Peru and Bolivia countries. Hence, researchers have been trying to predict the water level of the TL in accurate and cost-effective ways.
One of the main disadvantages of machine learning models is the difficulty in tuning model parameters for the optimal learning process since it influences the prediction performance of the models (Chen et al. 2017). Machine learning models are well known as hyper-parameter methodologies; hence, the optimization process for tuning their internal parameters is highly essential Mohammadi and Mehdizadeh 2020). One of the popular optimization procedures used in the literature is the trial and error experiment which is time-consuming, prone to data errors, and sometimes achieve unrealistic prediction (Leuenberger et al. 2018). In recent years, several studies have been conducted on the machine learning methods for predicting LWL in different regions (Coulibaly 2010;Buyukyildiz et al. 2014;Tongal and Berndtsson 2014;Aytek et al. 2014;Sanikhani et al. 2015;Kisi et al. 2015;Seo et al. 2015;Shiri et al. 2016;Vaheddoost et al. 2016;Shafaei and Kisi 2016;Yadav and Eliza 2017;Esbati et al. 2018;Ebtehaj et al. 2019;Ghorbani et al. 2019;Zaji et al. 2019;Wen et al. 2019;Bonakdari et al. 2019;Hrnjica and Bonacci 2019;Zhu et al. 2020;Yaseen et al. 2020). Tongal and Berndtsson (2014) applied a nonlinear approach, namely phase-space reconstruction and selfexciting threshold autoregressive model (SETAR), for forecasting daily measured LWL at three Sweden lakes, namely Vänern, Vättern, and Mälaren. The authors reported excellent results with a coefficient of correlation of 0.99 for the validation dataset. Aytek et al. (2014) investigated the capabilities and robustness of the gene expression programming (GEP) model for modeling monthly LWL using data measured at Van Lake, Turkey. The obtained results using the GEP model were compared to those obtained using the adaptive neuro-fuzzy inference system (ANFIS), and it was demonstrated that GEP was more accurate and achieved a coefficient of determination (R 2 ) of 0.989 compared to 0.972 obtained using the ANFIS model. Sanikhani et al. (2015) compared three data-driven models in forecasting monthly LWL of Manyas and Tuz lakes, Turkey. The proposed models were: (i) GEP, (ii) ANFIS with grid partition (ANFIS-G), and (iii) ANFIS with subtractive clustering (ANFIS-S). The three models were used for forecasting LWL at 1, 2, and 3 months in advance, and it was found that ANFIS-G was more accurate for forecasting LWL 1 month in advance with root mean square error (RMSE) of 0.251 and R 2 of 0.872, while ANFIS-S was more accurate for forecasting LWL at 2 and 3 months in advance. Esbati et al. (2018) compared the ANFIS model trained using four different training algorithms, namely: (i) ANFIS trained using the particle swarm optimization (ANFIS-PSO), (ii) ANFIS trained using the particle swarm optimization and recursive least square (ANFIS-PSORLS), (iii) ANFIS trained using hybrid backpropagation-RLS (ANFIS-BPRLS), and (iv) ANFIS trained using hybrid PSO and Gradient descent (ANFIS-PSOGD). The proposed models were applied and compared using the Urmia Lake, Iran, and the ANFIS-PSO data for forecasting LWL at 10-day ahead with R 2 of 0.986. A new bio-inspired algorithm called King's Castle Optimization (KCO) with training sample adaption (TSA) was proposed by Zaji et al. (2019) for optimizing the ANFIS parameters, and a new model called KCO-TSA-ANFIS was developed for forecasting monthly LWL of Lake Van, in Turkey. Results obtained show that the KCO-TSA-ANFIS was more accurate with R 2 and RMSE of 0.989 and 0.0602, compared to 0.973 and 0.208 obtained using the standard ANFIS. Hrnjica and Bonacci (2019) used the long short-term memory (LSTM) recurrent neural network) and the feedforward neural network (FFNN) for predicting the LWL of Vrana Lake in Croatia. The authors reported that LSTM was more accurate with a correlation coefficient of 0.897 compared to 0.852 obtained using the FFNN models. Wen et al. (2019) proposed applying the deep learning (DL) approach for predicting the LWL of Sumu Barun Jaran Lake, China.
The model was developed using several climatic variables, namely relative humidity, wind direction, precipitation, air temperature, and water temperature. Results obtained using the DL mode were compared to those obtained using the FFNN and the multiple linear regression (MLR), and the best accuracy was obtained using the DL (R 2 = 0.657), slightly higher than the FFNN (R 2 = 0.642) and significantly higher than the MLR model (R 2 = 0.554). Bonakdari et al. (2019) compared four machine learning models, namely extreme learning machine (ELM), Gaussian process regression (GPR), minimax probability machine regression (MPMR), and relevance vector machine (RVM), for predicting LWL for Michigan Huron lake located in the USA. By comparing six different input combinations, the authors reported that the ELM and MPMR were most accurate than RVM and GPR models. Yaseen et al. (2020) proposed a new hybrid bio-inspired model called whale optimization algorithm (WOA) optimized multilayer perceptron (MLP-WOA) for LWL forecasting. The hybrid new model was compared to standalone machine learning models, namely self-organizing map (SOM), FFNN, random forest (RF) regression, decision tree regression (DTR), and the cascade-correlation neural network model (CCNN). The models were applied and compared for forecasting LWL at Van Lake located in eastern Turkey, and the superiority of the MLP-WOA was particularly promising since they indicated the new approach superiority over the standalone machine learning. Finally, Zhu et al. (2020) used monthly data of 69 lakes in Poland for forecasting LWL one month in advance by applying two models: the DL and the FFNN. The authors argued that they found the DL slightly superior to the FFNN in only over one-half of the lakes, and its performances did not significantly surpass those of the FFNN model in terms of R and RMSE values. Ebtehaj et al. (2021) reported better performance of hybrid models (GS-GMDH: RMSE = 82.98 and ANFIS-FCM: RMSE = 83.8) compared to standalone models (GMDH: RMSE = 88.05 and GEP: RMSE = 88.22) for predicting water level at two stations in Malaysia.
The results of past research by Kisi et al. (2015) on Lake Urmia in the Middle East showed that bio-inspired algorithms are a reliable tool for estimating LWL. Their study coupled the firefly worm optimization algorithm with backbone regression coupling, which resulted in modeling with a coefficient of correlation of 0.99. Another study was conducted by Ghorbani et al. (2018) on Lake Egirdir in Turkey, using a hybrid firefly algorithm and an artificial neural network (ANN). The results of their study showed that the meta-search algorithms could increase the accuracy of the ANN model, and they also modeled the LWL of Lake Egirdir with an error of 0.148 m to 0.029 m. The water level of Lake Cypress and Lake Winnipesaukee in the USA was forecasted by Ghorbani et al. (2018). They have used a combined ANN approach and gravitational search algorithm and particle swarm optimization (PSO) techniques. They reported that the new hybrid models accurately forecasted LWL by RMSE ranging from 1.457 to 0.055 m.
Among several machine learning models, the ANFIS is one of the popular and remarkably adopted models for implementing hydrological process simulation (Ghose et al. 2013;Roushangar et al. 2014;Moazenzadeh et al. 2022;Hadadi et al. 2022). For this model, learning the membership functions and adjusting the connections between the layers has been the main challenge for the practitioners. The gradient descent approach is commonly used for ANFIS model parameters optimization. In this method, the gradient is calculated during several iterations. In each iteration, its performance was influenced by the initial point (Ewees and Elaziz 2018). In recent studies, metaheuristics optimization algorithms have proven their potential to efficiently overcome the difficulties encountered during the parameter tuning of machine learning models (Ewees and Elaziz 2018;Mohammadiet al. 2020). Metaheuristics optimization algorithms can help in estimating the parameters of the machine learning models automatically (Termeh et al. 2018) and improve the model's performance significantly. Various metaheuristics algorithms inspired by biological systems have been introduced recently to provide an efficient optimization approach by imitating various biological phenomena for hydrological goals. Some of those prevalent nature-inspired metaheuristics algorithms in hydrological forecasting include the genetic algorithm (Jahani and Mohammadi  Titicaca Lake is the most important and largest lake in South America. In addition, this lake plays a vital role in agriculture, the habitat for animals, human beings, and plant species and the ecology of the region, economic, industrial activities, and environment (Kroll et al. 2012). Titicaca Lake was designated as a Ramsar Site in August 1998 since it is home to several water birds (Pillco et al. 2019). The Titicaca Lake is an international wetland, the largest, deep, and the highest navigable water body globally; it is also very important from the point of view of the importance of tourism in South America. It is situated between Peru and Bolivia (between the latitude of 15°25' and 16°35' south). Titicaca Lake is the only home to multiple threatened species such as the huge Titicaca water frog and the flightless Titicaca grebe, and also a crucial surface pond for ecology in the western region of South America (Guédron et al. 2018). Titicaca Lake water level is a key component of the hydrological cycle in the Peruvian and Bolivian Altiplano. For planning, designing constructing, and operating lakeshore structures, the knowledge of LWL fluctuations is very crucial.
Furthermore, for the planning and management of freshwater supply, it is important to have accurate predictions of future LWL. Therefore, it is necessary to develop a simulation model of extreme or abnormal level variations in order to control future LWL changes. Accurate estimation of Titicaca LWL is a fundamental prerequisite for the sustainable management of water resources in South America and especially in Bolivia and Peru (Pillco et al. 2019).
This study is an important advancement for water resource planning and management in the region of Titicaca Lake in Bolivia and Peru. Firstly, the region of Titicaca Lake basin has experienced extreme climate change. Therefore, identifying surface water changes and behaviors of water reservoirs is a significant step in developing water resource management strategies. Secondly, the countries of Peru and Bolivia have become significant trading centers in coffee products in the world. Therefore, agricultural and environmental managers need to have precise planning in the water resources basin for agricultural activities. Planning water resources by using modern irrigation methods can help Peru and Bolivia to promote traditional agriculture to modern agriculture. Thirdly, this research will be useful with the newest methods of artificial intelligence combined with optimization methods to predict the highest accuracy of water resources, including the LWL. Fourthly, in this research, we have been focused on water resources in South America, and this research will help other researchers to get acquainted with the pattern of water resources in this particular region. With the motivation to allocate robust and reliable artificial intelligence models for hydrological process simulation, the current research emphasizes the implementation of the newly explored hybrid ANFIS model with a couple of metaheuristics algorithms. Three different optimization algorithms based on nature-inspired phenomena (i.e., GA, PSO, and WOA) are hybridized with the ANFIS model to improve predictability for monthly scale LWLs. The proposed hybrid models have been tested at Titicaca Lake for predicting monthly LWL, and their performance has been validated against basic seasonality evolutionary, the effect of monthly, and yearly lags on the models. The manuscript is structured as follows: Sect. 2 presents the study area, data description, and a brief explanation of the optimization algorithms and the predictive models. The experiment and results are discussed in Sect. 3. Finally, the conclusion and remarks are reported in the last section.

Study site and dataset explanation
Titicaca Lake (15°45' S, 69°25' W), also known as 'the highest navigable lake' globally, by volume of water and surface area, is the largest inland freshwater body in South America. This lake is located on the border of Peru and Bolivia in the northern parts of the endorheic Altiplano basin high in the Andes. The eastern and western sides of this lake are located in Bolivia and Peru, respectively. Titicaca Lake has a surface area of 8560 km 2 , maximum length, maximum width, maximum depth, average depth, and water volume of 190 km, 80 km, 282 m, 108 m, and 894 km 3 , respectively. Titicaca Lake (Roche et al. 1992) is fed by five major river systems (Ramis, Coata, Ilave, Huancané, and Suchez) in South America (Grove et al. 2003). There are other smaller streams (approximately 20), which also drain into Titicaca Lake. The 41 populated islands are situated in Titicaca Lake. Figure 1 shows the location of Titicaca Lake in the study area.
The values of LWL lag times data during the period of 1973-2017 were used to estimate LWL for the next months by the ANFIS, ANFIS-GA, ANFIS-PSO, and ANFIS-WOA models. Table 1 lists descriptive statistics of lake level data used. Total data was divided into two parts for training (calibration) and for testing the developed models. The first 75% of the data (from August 1974 to January 2007, i.e., a total of 390 records) was used for calibration, and the remaining 25% (from February 2007 to January 2017, i.e., a total of 120 records) was used for testing purpose. Monthly patterns of LWL for a representative Titicaca Lake are shown in Fig. 2. To better understand Titicaca LWL fluctuations conditions, we used separate monthly data series from some years during the study period (Fig. 3).

Adaptive neuro-fuzzy inference system (ANFIS)
ANFIS is a type of artificial neural network based on the Takagi-Sugeno fuzzy system. This method was developed in 1993 by Jang (1993). Since it integrates neural networks and fuzzy logic concepts, it can take advantage of both of them in a single frame (Ehteram et al. 2019). The appropriate structure of the ANFIS technique is selected according to the input data, membership degree, input and output membership rules and functions (Larrea et al. 2021). In the training phase, by modifying membership degree parameters based on acceptable error rates, input values are closer to real values. The ANFIS technique uses neural network learning algorithms and fuzzy logic to design nonlinear mapping between input and output space and has good capability in training, construction, and classification. It also has the advantage of allowing fuzzy laws to be extracted from numerical information or expert knowledge and to formally formulate a rule. In addition, it can regulate the complex transformation of human intelligence into fuzzy systems. Its inference system is in accordance with the set of fuzzy if-then laws that can be learned to approximate nonlinear functions. Therefore, ANFIS is considered a popular estimator tool in engineering fields. In this study, Fuzzy Cluster Means (FCM) method was used for fuzzy clustering based on the k-means clustering algorithm and the backpropagation method was used for finding the best optimization method in classical ANFIS ).

Genetic algorithm (GA)
The genetic algorithm (GA) is an iteration-based optimization and search technique based on genetics and natural selection rules. According to Holland, this algorithm performs the evolution on the gene symbol corresponding to the answers to a problem, just as the process of gene evolution takes place in nature (Jahani and Mohammadi 2019). This method starts with an initial guess, and in different periods the best answers, according to their performance in the objective function and meeting the criteria and constraints, are selected and transferred to the next period. This algorithm considers several points of the search space simultaneously, thus reducing the chances of converging to a local minimum (a risk that limits most conventional search methods). Because in this algorithm, the value of each point is calculated independently, the simultaneous calculation of several points in the answer space can be performed using parallel computers. This algorithm does not require much  information about the problem and even its derivative information because this method in each period by calculating the value of each variable and generates new sets using its operators to generate sets of variables that create the best value in the problem and by repeating this action, it finally reaches a set of optimal points. GA uses several functions, which are described below. Selection operator: This operator selects a number of chromosomes from among the chromosomes in a population to reproduce. More graceful chromosomes are more likely to be selected for reproduction. Merge Operator: During the merger operation, parts of the chromosomes are randomly exchanged. Mutation operator: After the integration, the operation is completed, and the mutation operator acts on the chromosomes (Jahani and Mohammadi 2019). This operator randomly selects a gene from a chromosome and then changes the content of that gene. If the gene is of the binary number genus, it converts it to inversion, and if it belongs to a set, it replaces that gene with another value or element of that set.

Particle swarm optimization (PSO)
PSO metaheuristic algorithm is firstly proposed by Kennedy and Eberhart (1995), and it is designed to solve nonlinear optimization problems with continuous variables. In addition, unlike other evolutionary methods such as genetic algorithms, PSO can be implemented with a more straightforward program. This capability of PSO is one of its advantages compared to other optimization methods. The principles of the particle clustering algorithms are based on the fact that individuals in each group use two types of information in the decision-making process: first, their own experience, and second, the experiences of others (Assareh et al. 2010).
The PSO algorithm starts with a set of random answers (Swarm). Each member of this set is called a particle. Particle navigation is done in such a way that the entire particles store in their memory the best position they have gained during the search process called Pbest. On the other hand, the best position achieved by each particle up to each stage is maintained under the name of Gbest. In this whole algorithm which is based on a weighted average of random components, the particles move toward better solutions, such as Pbest and Gbest, to eventually converge to a single point. Position (X) and velocity (V) are the main parameters of this algorithm which can be explained as follows ): where w is the coefficient of inertia, r 1 and r 2 are random coefficients, c 1 and c 2 are the acceleration coefficients, X pbesti and X gbesti are the best position of ith particle and global value in each iteration, respectively. V i (t + 1) and X i (t + 1) are the next velocity and new position, respectively.

Whale optimization algorithm (WOA)
The WOA as a meta-metaheuristic optimization algorithm is one of the newest population-based optimization algorithms, which was first proposed in 2016 by Mirjalili and Lewis (Mirjalili and Lewis 2016). The WOA is inspired by the nature and social behavior of whales; this algorithm uses a bubble net hunting strategy to explore and exploit, and by avoiding local optimal points, with an integrated adaptive technique, can achieve the optimal solution with less computational time wasted (Soleimanian Gharehchopogh and Gholizadeh 2019; Vaheddoost et al. 2020). The WOA is performed in three stages or three phases, which are as follows: (1) Siege hunting; (2) Operation phase: method of attacking the bubble net; (3) Exploration stage: Hunting search. The WOA algorithm starts with a set of random solutions. The search agents update their position in each iteration according to the randomly selected search agent with the best current solution (Mohammadi and Mehdizadeh 2020). (1)

Modeling development
The goal of this study is to examine the feasibility and to demonstrate the benefits of improving the basic ANFIS structure using metaheuristics algorithms (i.e., ANFIS-GA, ANFIS-PSO, and ANFIS-WOA) for providing a rapid and efficient LWL prediction method at the Titicaca Lake. The performances of the GA, PSO, and WOA algorithms were examined, and the weight and bias of the basic ANFIS model were optimized by the bio-inspired optimization approaches. The hybrid models (i.e., ANFIS-GA, ANFIS-PSO, and ANFIS-WOA) stop when a mathematical fit between ANFIS weights and the GA, PSO, and WOA is created or the maximum number of iterations occurs. It is an estimator hybrid procedure that utilizes both capabilities of the ANFIS and optimization algorithms. Some research has shown that such a hybrid technique can present more successful results (Ghorbani et al. 2018;Mohammadi and Mehdizadeh 2020;Mohammadi and Aghashariatmadari 2020;Pham et al. 2020).
The proposed modeling method is shown in Fig. 4, starting by normalizing the dataset as in Eq. (6), then after, lag months based on autocorrelations were employed to identify the most robust and efficient model for monthly scale LWL simulations. A scenario with the most correlation value is assigned as the most suitable input of models for LWL simulation. Three hybrid models are used for this aim; single ANFIS trained by GA, PSO, and WOA. As one of the advantages of using a hybrid model is these three optimization techniques improved the learning system of classical ANFIS, and then it gets better performance.

Data normalization
Data is normalized by scaling the various time series related to the input and output variables within a specified range. In this research, the data time series were normalized in the range 0 and + 1 by using the following equation: where X n explains the normalized input values lying in the range of [− 1, + 1], X i is the observed input values, and X min and X max explain the minimum and maximum values of the explanatory variables.

Model evaluation
Various performance evaluation criteria such as coefficient of determination (R-squared) (Essam et al. 2022), mean absolute error (MAE) (Ahmed et al. 2018), and root mean square error (RMSE) (Abdullah et al. 2019) were utilized to evaluate the simulated results of each model, as illustrated in Eqs. (4), (5), and (6), respectively: where n represents the total length of the evaluated data, x and y, respectively, demonstrate the observed and simulated data and σ x and σ y are the standard deviations of the observed and simulated data.

Model identification
The predictive accuracy of the machine learning approaches is significantly related to the proper selection of the input combination and efficient size of the training data (Aghelpour et al. 2019). In this regard, the autocorrelation technique is used to lag months and the effect of seasonality on LWL for finding the best mask of the input variables. The diagram of the Titicaca LWL autocorrelation from 1973 to n 2017 is shown in Fig. 5. As shown in Fig. 5, the water level of Titicaca Lake is affected by seasons.
In other words, seasonality has a significant impact on the water level of the lake. Therefore, the dynamic seasonal memory of the LWL has been used as a predictor model. So, four scenarios with short-term delays (1, 2, 3, and 4 months) as well as seasonal long-term (12, 24, 36, and 48 months) have been defined, and these scenarios are listed in Table 2. The variable combinations for each of the defined scenarios are used as inputs for ANFIS, ANFIS-GA, ANFIS-PSO, and ANFIS-WOA models to predict the LWL of Titicaca Lake.

Result of models' error estimation
The ANFIS and hybrid ANFIS models optimized using metaheuristic algorithms (including GA, PSO, and WO) have been implemented with respect to the predetermined inputs, and the results are shown in Table 3. The best ANFIS output was obtained using the ANFIS2 scenario, where the RMSE, MAE, and R 2 values were 0.136 m, 0.108 m, and 0.925, respectively; hence, the ANFIS2 has been chosen as the best model for ANFIS base models (ANFIS1, ANFIS2, ANFIS3, and ANFIS4). However, in ANFIS hybrid models, changing the input scenario causes very little change in model accuracy. Therefore, in such situations where the forecast error difference is about 0.001 m, the first scenario (ANFIS-GA1, ANFIS-PSO1, and ANFIS-WOA1), which includes the lowest number of input variables, can be selected for best (high accuracy) models. The reason for this can be the complexity of the model's structure, which means that the complexity of the hybrid models makes it possible to achieve the desired result with the lowest input. However, in the ordinary ANFIS model, in order to increase the accuracy of prediction algorithms, we should pay particular attention to the increase of input variables.
Comparison of the results among the models shows that GA, PSO, and WOA optimization algorithms have enhanced the prediction accuracy of the ordinary ANFIS model as follows: 57, 54, and 70 percent (%), respectively, in the first scenario, 41, 43, and 51 percent (%), respectively, in the second scenario, 47, 47, and 55 percent (%), respectively, in the third scenario, and in the fourth scenario 53, 53, and 55 percent (%). The worst performing algorithm was PSO, followed by GA and WOA, respectively. ANFIS-WOA1 model, which had the best prediction accuracy, was found to have an RMSE value of 0.087 m, an MAE value of 0.063 m, and R 2 of 0.964.
The correlation of the best predictions of each model compared with the observed data is shown for the training and testing phases in Fig. 6. It is clearly evident from Fig. 6 that the outputs of each of the four models are assembled in a satisfactory way around the 1:1 line, but generally, the distance from the 1:1 line is larger in the ANFIS outputs. The points in the scatter plots are closer to the 1:1 line for ANFIS-GA and ANFIS-PSO models, which implies that the calculated values are closely matching with the observational values for these models. The best precision around the 1:1 line is also for the prediction of ANFIS-WOA, which means that even some relatively remote points in this model have reached their nearest distance to the 1:1 line. Another point in these graphs that needs to be highlighted is the higher correlation between outputs during the training period than the test period, which indicates satisfactory model training (calibration). As well-trained models can have better estimates of the lake's water level, and ultimately, with better training models, it is possible to obtain higher accuracy in estimating the LWL. Also, the range of changes in the LWL during the test period is significantly smaller than the training period, so the reported error in the test period also has a smaller range of amplitudes.
In this section, the violin plot has been used to compare the distribution of the observed and predicted LWL datasets in the testing period (Fig. 7). The violin plot has shown that the frequency distribution of the outputs of hybrid models is not significantly different at different Lake level(X) t ANFIS2, ANFIS-GA2, ANFIS-PSO2, ANFIS-WOA2 Lake level(X) t ANFIS4, ANFIS-GA4, ANFIS-PSO4, ANFIS-WOA4 inputs, and the same distribution is reported for prediction provided by all four input patterns. Therefore, due to lower inputs and the estimation of optimal frequency distribution, input model number one can be considered as superior for hybrid models. Different inputs for ANFIS are somewhat impressive; for example, the ANFIS model has the weakest median score in the first scenario (by inputs of X t−12 ). The ANFIS model in the second scenario (by inputs of X t−1 , X t−2 , X t−3 , X t−4 , X t−12 , X t−24 ) estimated the first quartile at about 3808.21 m (first quartile of observation values = 3808.31 m), in the third scenario (by inputs of X t−1 , X t−2 , X t−3 , X t−4 , X t−12 , X t−24 , X t−36 ) estimated the third quartile about 3809.06 m (third quartile of observation values = 3808.98 m), which is the weakest function in estimating quartiles, among all the predictions. In hybrid models, the median and quartile estimates are very close to real values. The estimate of the frequency distribution of data has improved in the ANFIS-GA and ANFIS-PSO models to a large extent. However, in the ANFIS-WOA model, the curvature of the violin is very close to the curvature of the violin of observational values. In Fig. 8, the best output of each model has been plotted along with the measured values of the time. This chart is drawn up to understand the ability of each model better for the test period. Also, Fig. 8 shows an acceptable overlap in the predictions provided by the applied models. The ANFIS model had underestimated in some spots, including March-August 2013, January-May 2014, August-October 2014, January-February 2015, and August-October 2015. The hybrid models have better overlap with observational values, as their disagreements with observational values are significantly less than the ordinary ANFIS model. The only disadvantage of hybrid models is that they had overestimated in April months; for clarity, in Fig. 8, the 2 years 2008 and 2016 have been isolated and exhibited. Most of the overestimated errors were in the April estimates for ANFIS-GA and ANFIS-PSO outputs, and it reported the lowest for ANFIS-WOA. Also, the models' parameters are listed in Table 4.
Further, Zhu et al. (2020) suggested that machine learning models are generally good at predicting the LWL; however, special care should be taken at the calibrating stage of the models. Hrrnjica and Bonacci (2019) suggested that the ANN models generally outperform the traditional ARIMA models in predicting LWL. These reported results validate our obtained results and prove the superiority of the machine learning models over traditional techniques. Choi et al. (2020) assessed the various machine learning models and found that more advanced models are superior as compared to traditional ANN models. In recent years, machine learning approaches have emerged as an alternative to complex physical-based models that need significant data and computation resources for calibration. Therefore, there is a need to test different varieties of these algorithms to test and establish the reliability of these methods for various hydrological applications. This study tested several metaheuristics optimization-based machine learning algorithms to test their accuracy in LWL prediction. The results of this study suggest that all methods performed very well for LWL prediction and can be utilized for field-level prediction. However, the WOA method was found to be better than the others for LWL prediction in our study. A similar model structure can be used for other lakes in different parts of the world. However, rigorous testing should be done to verify the reliability of the adapted method.
In the literature, we can find several studies where similar algorithms were checked for various hydrological applications and different results were obtained. As per the WOA technique's ability for optimizing results of predictor tools such as ANFIS, MLP, and SVR, these results are similar to the result of Mohammadi and Mehdizadeh (2020), Vaheddoost et al. (2020), and Diop et al. (2020) as they reported ability of WOA techniques in improving the accuracy of hydrological variables modeling. Diop et al (2020) used WOA combined with ANN to predict rainfall and found that WOA optimization-based results are better than other models and WOA improves the prediction accuracy. Similar results were also obtained by Vaheddoost et al. (2020) for estimating field capacity and permanent wilting points of the soil. These studies clearly confirm our findings since, in our study; also we have found that the ANFIS-WOA was the best approach. However, in this study, we demonstrated that the proposed hybrid model results in a better estimation of LWL. Multiple previous studies have reported that, in general, hybrid models perform better than simpler models due to their increased complexity and nonlinearity (Chen 2013;Pham et al. 2019;Mohammadi and Mehdizadeh 2020). However, with an increase in complexity and nonlinearity, significant care should be taken to avoid overfitting (Kim and Park 2006;Vaheddoost et al. 2020).
The results of this study indicate that complex and advanced machine learning models based on metaheuristics approaches are having significant potential and are generally better than the traditional machine learning models. The

Conclusions and recommendations
Titicaca Lake is one of the most important water bodies for water supply in the South American continent. Better forecasting algorithms for LWL fluctuation will substantially help in operational management and better development of infrastructure in the region. It has a total surface of 8560 km 2 , and it is the largest and deep navigable water body in the world. So, this study proposed new meta-optimized hybrid models to simulate the monthly water levels of the lake. The lag times of the Titicaca Lake water level from Bolivia and Peru were used to implement the models. The swarm-based GA, PSO, and WOA metaheuristic optimization algorithms were utilized for optimizing the ANFIS model parameters, and the results were compared with those of the basic ANFIS model. The results demonstrated that meta-optimized hybrid models enhanced the capability of the original ANFIS model for the reproduction of monthly scale LWL. The meta-optimized model also outperformed the evolutionary-based Takagi-Sugeno function. The ANFIS is the predictive model for this research, and the Takagi-Sugeno function was employed to learn the original ANFIS model parameters. Among the meta-optimized models, the ANFIS-WOA model was found to be the most suitable model because of the lowest total error rate and fastest convergence speed. The ANFIS-WOA model still showed comparable results to ANFIS-PSO in accurately reproducing the high flows. The WOA algorithm has superiority over the other metaheuristic algorithms because of its capability of escaping from the local optima and its faster convergence speed. The capability of this algorithm to escape from local optima is due to the harmony between the search and exploitation phases, which distinguishes it from other metaheuristic algorithms. The results indicate that the capability of WOA for enhancing the forecasting accuracy of lake water-level models. Furthermore, this work highlighted the superiority of hybrid models over other simple models that can be extended for other engineering applications in developing the predictive model for various features. Finally, the ANFIS-WOA was utilized for monthly scale LWL simulation in Titicaca Lake. The calibration and verification of ANFIS-WOA were performed against the historical monthly LWL data. The statistical results of the performance evaluation of the proposed models are also shown, when used to learn the ANFIS model parameters, the WOA algorithm improved the ANFIS model efficiency in terms of RMSE, MAE, and R 2 statistical measures. The ANFIS-WOA model for the best scenario (inputs were X t−1 , X t−2 , X t−3 , X t−4 , X t−12 ) with the RMSE = 0.087 m, MAE = 0.063 m, and R 2 = 0.96, outperformed the original ANFIS model for the best scenario (inputs were X t−1 , X t−2 , X t−3 , X t−4 , X t−12 , X t−24 ) with RMSE ≈ 0.15 m, MAE ≈ 0.11 m, and R 2 ≈ 0.91. The WOA trainer also outperformed other metaheuristic trainers, such as GA and PSO algorithms. The model inherits the high exploration capabilities of the WOA algorithms to find the best solution. This study indicated the robustness of the WOA-based trainer to simulate the monthly scale LWL data. However, further studies are encouraged to evaluate this model in different hydrology regions with diverse complexity levels of hydrological data. A higher temporal resolution model (e.g., daily or hourly discharge) can also be tested in the future to extend this study, which is likely to generate more robust models for higher resolution prediction. otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.