Monitoring network design with MLPG-TLBO hybrid model (case study Birjand, Iran)

As the groundwater quantitative monitoring aimed to determine the factors affecting the aquifer behavior plays an important role in its regional management, studying the temporal and spatial groundwater level variations requires a comprehensive monitoring network. Effort has been made in this study to introduce a new linked simulation–optimization method, named MLPG-TLBO to quantitatively monitor the Birjand aquifer and determine the optimal points for piezometers. This model uses meshless local Petrov Galerkin (MLPF) method in the simulation part and teaching–learning-based optimization (TLBO) method in the optimization part. The objective function, in this study, is to minimize the difference between the groundwater level observed in piezometers and obtained computationally by the model. Since this coupled model is independent from the meshing, it eliminates the mesh-dependent shortcomings and, hence, yields more accurate results. It has been calibrated and validated in previous Birjand area studies and has led to acceptable results. By implementing the model in Birjand aquifer, the optimal positions of ten piezometers were determined mostly in areas where the density of the extraction wells was lower. Finally, the RMSE of both MLPG-TLBO and FDM was obtained to be 0.334 m and 1.483 m for 10 optimal piezometers. The RMSE value for MLPG-TLBO has shown a 0.423 m reduction compared to its previous value. This difference is quite meaningful as it shows good performance of this method in designing an optimal network for the aquifer.


Introduction
The hydrometric monitoring network plays an important role in the planning and management of water resources and its design for the groundwater aquifer can either be qualitative or quantitative, but this does not mean that they do not do each other's duty; it rather shows that their design methods change according to the users' objectives. Therefore, in a monitoring network proportionate to the aquifer conditions, the concepts of protection and maintenance, as well as the development and updating to receive highly uncertain hydrometric information, will be very important (Mishara and Coulibaly 2009).
The hydrogeological modeling methods of the groundwater monitoring are either statistical or probabilistic, both of which require the access to sufficient information on the quantitative and qualitative status of the aquifer for optimal management (Mohtashami et al. 2017). The groundwater level study/evaluation through the existing models facilitates the behavior prediction of groundwater resources. Since the groundwater monitoring network design is an optimization issue requiring high computation costs, it should provide the best results with the least number of observation wells.
Studies on the quantitative and qualitative design of the aquifer monitoring network are numerous. Using the genetic algorithm, Cieniawski et al. (1995) determined the location of groundwater monitoring wells to examine such objectives as the maximum reliability and minimum contaminated areas in the initial diagnoses. They considered 162 monitoring-well positions, evaluated the Pareto front values for 4 and 5, and showed that the reliability was 48-90% and contamination was 3.5-14.5 hectares based on the considered conditions and results (Cieniawski et al. 1995). In a study performed in the California National Laboratory for groundwater quality monitoring, Kollat et al. (2011) used the NSGAII algorithm aimed to reduce the contaminant concentration estimation errors, contamination load estimation errors and costs caused by the leakage of the PCE contaminant from an underground tank (Kollat et al. 2011). There were 29 wells in the flow path, and the Pareto front was obtained for 18 to 25 monitoring wells. They showed that the algorithm was able, in all the considered cases, to provide different options for the location of monitoring wells in accordance with the objective functions. To optimally design the groundwater quality monitoring network and show spatiotemporal densities, Dhar and Patil (2012) used the non-dominated sorting multi-objective genetic optimization algorithm and the conventional fuzzy Kriging method to determine the minimum number of piezometric wells in a real Aquifer of India (Dhar and Patil 2012). Aboutalebi et al. (2017) used the NSGAII algorithm for the optimal design of monitoring networks in river-reservoir systems and minimized the average contaminant detection time and concentration prediction error at sampling stations using the multi-objective optimization method (Aboutalebi et al. 2017). Aiming at reducing the number of piezometers without affecting the data quality and using KGA (Kriging-combined genetic algorithm) and PSO (particle swarm optimization) algorithms, Khorramdel and Kekhaei (2016) studied the Astaneh-Kuchesfahan aquifer in the north of Iran to optimize the monitoring network. They showed, through the results, that the number of observation wells in this aquifer could reduce by 26% without any losses in the data quality (Khorramdel and Kekhaei 2016). Aiming at minimizing the number of piezometric wells and estimation error, Mirzaei Nodooshan et al. (2016) used the NSGAII algorithm to locate piezometric wells and introduced 311 points in Eshtehard aquifer to the NSGAII algorithm as potential points for which the observational groundwater levels were interpolated using real piezometers and interpolation tools in the GIS environment. Locating was done by the NSGAII algorithm and the decision variable was the piezometric wells (Mirzaei Nodooshan et al. 2016). Mirzaei Nodooshan et al. (2019 presented an optimization-based method for the network design of the groundwater level monitoring and used it for the Eshteharud plain aquifer. Scenarios they considered to optimize the position of observational wells were: (1) Selecting optimum number of wells from among those in the study area and (2) Adding optimum number of monitoring wells to the existing ones. Results showed that their method could omit a number of extra wells in an existing monitoring network with acceptable accuracy or add some if it lacked any (Mirzaei Nodooshan et al. 2019). To monitor the Gilan plain aquifer qualitatively, Janatrostami and Salahi (2020) used the GA optimization method and showed that it was acceptable based on the calculated error criteria. Since the spatial distribution studies and comparison of the average observed regional salinity with that calculated in the optimal network (covering all selected monitoring wells) showed that the calculated values were larger than the observed average, the network provided the optimal data of the groundwater quality in more polluted areas (Janatrostami and Salahi 2020). Temur (2021) used the hybrid TLBO algorithm for Optimum design of cantilever retaining walls under seismic loads. The robust results were obtained by HTLBO with 50 populations. Consequently, the results showed that the hybrid TLBO algorithm optimized the implementation cost of RC cantilever retaining walls well (Temur 2021). Recently, various researchers have used combined optimization algorithms to overcome computational complexity. Veladi and Beig Zali (2021), to optimally design single-layer dome structures, used the hybrid TLBO algorithm. The analysis results during this comparative process indicate that the hybrid TLBO algorithm provides a better and simpler design in less time than the standard algorithms (Veladi and Beig Zali 2021).
A literature review reveals that research that combines numerical and optimization models to design a monitoring network is rare. Numerical models simulate the desired phenomenon with great accuracy and predict and simulate its behavior through the provided input information (Khorashadizadeh et al. 2020). Those used in the groundwater domain are summarized in the finite difference and finite element methods, which form the basis for GMS and FEFLOW, respectively (Majidi Khalilabad et al. 2021). Use is recently made, in this domain, of new meshless numerical methods that can be improved by the user-required possibilities besides being more accurate and capable than their predecessors due to their being open-source (Liu and Gu 2001). Meshless methods have a wide range including the "smoothed-particle hydrodynamics (SPH)," "elementfree Galerkin (EFG)," "meshless local Petrov-Galerkin (MLPG)," and so on. As choosing a method depends on its application, it is better to use specific methods in fluidrelated problems.
The aim of this study is to use a numerical and an optimization model together to design a monitoring network optimally. The MLPG numerical method and the TLBO (training-learning-based optimization) algorithm were linked together recognized as simulation and optimization models to develop a simulation-optimization model for the design of an optimal monitoring network on Birjand aquifer in South Khorasan province.

Case study
Birjand aquifer, in South Khorasan province, is 269 km 2 in area containing 190 extraction wells and 10 piezometric wells. It is 1491 m above the sea level with an annual rainfall of 160 mm and is known as a dry-climate region according to the Domarton classification. Figure 1 shows the location of the aquifer, extraction wells (blue) and piezometric wells (red) (Khorashadizadeh et al. 2016 and. Rao et al. (2011) were the first to introduce the TLBO algorithm in the "Computer-Aided Design" Journal to solve the related problems (Rao et al. 2011), but about a year later, they introduced it in the "Information Sciences" Journal in a larger scale and for more general problems (Rao et al. 2012). Contrary to other meta-heuristic methods that usually model the nature, the TLBO algorithm (described later in more details) models the classroom teacher-student process and uses it to find the optimal solution. It is based on teaching (teacher) and learning (student) steps and its output at each step is the students' scores and level of knowledge (Črepinšek et al. 2012).

Generation of initial population
In the first step, a series of points (students) are formed and distributed randomly in the solution space and then, the objective function (students grades and level of knowledge) is calculated for each point based on the problem's main objective function and coordinates of each point; the score distribution has been assumed to be bell-shaped (Fig. 2) without reducing the problem all-inclusiveness. As shown, the best student becomes the teacher and the training step begins.

Teaching stage
Here, the point specified the "best" by the objective function becomes the teacher and starts teaching others; in other words, all the points with situations worse than the teacher move toward the latter causing the average distribution to move toward it. It is worth noting that in meta-heuristic algorithms, since use is always made of a randomizing factor to randomly check the points, the distribution movement toward the teacher, too, has the mentioned factor. In short, for any distribution point x i , there is a relation as follows: where r i is a random number between 0 and 1, T A is the coordinates of the teacher point, and T f is the teacher training factor selected randomly to be either 1 or 2. (1) has a higher value than x i , it will replace it. Training step changes are shown in Fig. 3.

Learning stage
Next is the learning step where the mutual learning of the students of one class is simulated in such a way that for each student, another one is randomly selected. Then, if the second has a better position than the first as regards the objective function, the first will move toward the second; otherwise, the second will move toward the first.
This algorithm-diversifying move too has a random factor, and the step for each student is as follows: Next is the checking of algorithm termination conditions; if conditions are met, the algorithm terminates, and if not, the loop is repeated.

Meshless local Petrov-Galerkin (MLPG)
Used for the first time by Atluri and Zhu (1998) to solve potential equations (Atluri and Zhu 1998), MLPG (meshless local Petrov-Galerkin) has since been used by a large number of researchers in scientific and engineering problems. In this method, since the test and experiment functions are selected from two different spaces, many meshless methods have used it as a basis for their formulation (Akbari et al. 2010). It uses the moving least squares approximation function to construct the shape function, but as the latter does not satisfy the Delta Kronecker function conditions, it fails to apply the problem's necessary boundary conditions (Dirichlet Boundary Conditions) making the use of such other methods as the "collocation" and "penalty function" a necessity (Liu and Gu 2005). But since these methods add approximate values to the results, it would be better to replace with those that eliminate these approximations and bring the model Fig. 2 Bell distribution of students grads Fig. 3 Changes in grads at teaching stage closer to reality. Although error rates are small in these methods, their removal and replacement with the moving Kriging method is a serious measure because by satisfying the Delta Kronecker function, it can apply the necessary boundary conditions without adding an approximation . As it is necessary in this research to study and analyze the uncertainty of the problem's boundary conditions, the moving Kriging method is used to directly apply the necessary ones. According to different references, this method has been used since 2013 to do approximation and construct shape function in the field of solids; however, this is the first time it is being used in this study in the field of fluids, especially groundwater (Khankham et al. 2015).

Moving kriging
The moving Kriging approximation function is a wellknown method for spatial interpolation in geometry and mining (Khankham et al. 2015 Matrices A and B are defined in Eqs. (4) and (5): where I is an N × N unit matrix and vector p (x) is as follows: In N × M matrix P, base polynomials are written for each point as follows: Matrix R and vector r(x) are found as follows: where r x i , x j is a correlation function between two points located at X i and X j and represents the covariance of H(x) real value: The above equations show that matrix R and r x i , x j play important roles in the calculations. The simple correlation function, used frequently, is a Gaussian function : where r ij = x i − x j and d c and c > 0 are constant-value correlation parameters used for fitting.

Governed equation of groundwater in unconfined aquifer
When the aquifer is unconfined, the governing equation in an unsteady state will be as follows based on Dupuit assumption and continuity equations (Dupuit 1863;Mays and Todd 2005): where H is the groundwater level (m), k is the hydraulic conductivity (m/day), s k x is in the horizontal direction, k y . is in the vertical direction, S y . is dimensionless specific discharge, Q is concentrated discharge (− ive if exploitation well and + ive, m 3 /day, if injection well), and q is the distributed flow (+ ive, e.g. rain or − ive, e.g. evaporation) (Mohtashami et al. 2020). It is worth noting that since the equation is time-depdent (unsteady state), the initial conditions are derived from the steady state and then put in the equation (Mohtashami et al. 2019).

Objective function (OF)
In this study, the objective functions are defined as in Eqs. (14) and (15) borrowed from (Jafarzadeh and Khashei-Siuki (2018), Cieniawski et al. (1995): In Eq. (14), h Obs . and h Sim are water levels observed in piezometers and simulated by the simulator model, respectively, E is estimation error, and F is the value of the objective function. The selected objective function varies in a 0-1 range and error E > 2.3 m will cause the objective function to tend to 1. Applying this objective function, especially for estimating the groundwater level, has caused the algorithm to search and discover the optimal points with very high severity.
The flowchart diagram of the MLPG-TLBO model is shown below (Fig. 4). The model, as mentioned in the introduction, has been formed by connecting two simulation and optimization models. The left part in the diagram shows the optimization section run athe main program of the simulation model. In the optimization diagram, the simulation model is executed after the parameters of the optimization model are given their initial values to calculate the objective function; then, the optimization steps continue after the initial values of the objective function are obtained (Table 1).

Results and discussion
To model the groundwater flow by the MLPG method, Birjand aquifer was first drawn by spreading points on its slope and boundary 500 m apart (Fig. 5) based on its earlier studies (Mohtashami et al. 2017) which enable the results obtained with this method to be compared with those of (Sadeghi tabas et al. 2016) and (Hamraz et al. 2015) based on the finite difference method. Next, the hydrodynamic information of each point was assigned to that node according to its position in the aquifer the hydrodynamic parameters of which include hydraulic conductivity and specific discharge information calibrated and validated by ( Sadeghi tabas et al. 2016). The aquifer boundary conditions defined for the related points (shown by arrows in Fig. 1) are 0 for constant head and flow. It is worth mentioning that this model was calibrated and validated by Mohtashami et al. (2017) in both steady and unsteady modes the results of which are listed in Table 2; modeling period was one year with monthly time-steps.
As optimization requires the piezometers' observational values to be specified throughout the aquifer and for each time-step, they were interpolated with the ArcGIS software in each month and at the aquifer surface, and an observational value was assigned to each point in each time period. Then, the objective function was calculated based on these values and those obtained from the simulator model and checked, based on the minimum value of this function, to determine the best location for the piezometric well in the aquifer. Table 2, the results of which were calibrated and validated based on the research of Mohtashami et al. (2017), compares errors of the MLPG model and those of the finite difference method.
The authors of this paper used the MLPG method proposed in previous studies to prepare, calibrate and validate the groundwater flow model of the Birjand aquifer. They had already calibrated the model in both steady and unsteady states and compared the results with those of the finite difference method (MODFLOW Package). The better and more accurate results of the meshless method were an evidence of the superiority of this model over the limited difference method due to its independence from the aquifer meshing which led to the elimination of the related errors. Results of these studies are presented in Table 2.
After coupling the meshless model to the optimization algorithm, the simulation-optimization calculations were done on a CORE I5/8 RAM/5200 U PC. This hybrid model was implemented 10 series the best of which was that of series 9 (with an objective function value = 0.132 m). The model implementation results are shown in Table 3.
Finally, the optimal location of the observational wells was selected based on series 9 implementation of the model. Piezometers' new positions are marked in blue in Fig. 6 where green indicates the previous locations of the observational wells in the aquifer.
As shown in the aquifer area, the optimal observational wells are located where extraction wells are less concentrated. While the spatial variations of wells 4, 6 and 9 are very high, those of 1, 2, 3, 5, 7 and 10 are almost as much as one node due to the proximity of observational wells to areas where the density of extraction wells is lower; this is quite checkable in well No. 4. New locations of observational wells are shown in Table 4 with their UTM coordinates.
The simulated and observed water levels in new observational wells are plotted on a graph and presented for piezometers 1, 4 and 7 in Fig. 7 to evaluate the performance of the optimization model to find the optimal location of the observational wells.
Water levels observed in new piezometers found by the meshless model and calculated by the finite difference method are compared in Fig. 7. In graphs of piezometers 1, 4, and 7 at their new locations in the aquifer, black, red and gray lines represent the water levels observed during modeling period, modeled by the meshless model and found by the finite difference method, respectively. In all diagrams, the meshless model yielded better and more accurate results than the finite difference method due to its independence from the meshing of the solution domain. In Fig. 7a (piezometer 1), the meshless model and observational results do not differ much at the beginning of the period, but over time, the difference increases gradually and reaches 0.07 m in the last month; for the finite difference model, this difference is  Fig. 7b, results of the meshless model are quite close to the observational data almost correlating until the month 6. Another point in 7b is the increase in the water level results of the finite difference method observable in periods 2 and 3 due to method-related errors. In Fig. 7c, unlike 7a, the meshless model is more accurate at the end of the period than at the beginning. Figure 8 shows the regression coefficient of the results of the meshless and finite difference models where green dots represent the water level obtained from the meshless model in new piezometers and the orange dots show those of the finite difference model. As shown, the regression coefficient of the meshless model is higher than that of the finite difference model in all piezometers; the higher is this coefficient, the greater is the correlation between the model results and the observational data.
Finally, to evaluate the model accuracy, several known error criteria in the groundwater field were used and the mean error (ME), mean absolute error (MAE) and root mean square error (RMSE) indices were found for all piezometers (Table 5).   As shown, displacing piezometric wells to the optimal positions found by the MLPG-TLBO model reduced the RMSE by 0.423 m which is quite high and meaningful for such an aquifer volume and means the monitoring network designed to operate in this aquifer is suitable.

Conclusions
This study proposed a novel model, i.e. MLPG-TLBO, as a simulation-optimization method aiming at enhancing the accuracy of groundwater monitoring and determining the optimal locations for piezometers in Birjand aquifer. After coupling the teaching-learning-based optimization model to the meshless flow model, already calibrated and validated in some studies, the optimal locations of ten piezometers were determined in Birjand aquifer and their observed water levels were compared with those found in different modeling periods computationally to evaluate the model performance proved to be acceptable because before applying the model, the RMSE was 0.757 m, but after optimization it reached 0.334 m. Considering the aquifer area, this difference is a huge volume of water the error of which was omitted by the MLPG-TLBO model. The acquired results also show the satisfactory performance of this model in this aquifer. It is noticeable that, in addition of giving a pervasive knowledge about groundwater behavior, this model helps the water engineers to define some scenarios for better management.   Funding The author(s) received no specific funding for this work not applicable.

Data Availability
The datasets used during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.

Consent for publication
This manuscript also does not contain data from any individual person, and therefore, consent to publish is not applicable.