Inverse Modeling of Seepage Parameters Based on an Improved Gray Wolf Optimizer

Abstract

The seepage parameters of the dam body and dam foundation are difficult to determine accurately and quickly. Based on the inverse analysis, a Gray Wolf Optimizer (GWO) was introduced into this study to search the target hydraulic conductivity. A novel approach for initialization, a polynomial-based nonlinear convergence factor, and weighting factors based on Euclidean norms and hierarchy were applied to improve GWO. The practicability and effectiveness of Improved Gray Wolf Optimizer (IGWO) were evaluated by numerical experiments. Taking Kakiwa dam located on the Muli River of China as a case, an inversion analysis for seepage parameters was accomplished by adopting the proposed optimization algorithm. The simulated hydraulic heads and seepage volume agree with measurements obtained from piezometers and measuring weir. The steady seepage field of the dam was analyzed. The results indicate the feasibility of IGWO in determining the seepage parameters of Kakiwa dam.

1. Introduction

In hydraulic engineering, seepage parameters of dams and dam foundations change with operating time and loading conditions. The changes in seepage parameters weaken the strength of the structure and lead to failure. Seepage analysis is commonly used to monitor the working conditions of dams and dam foundations for the safety of hydraulic projects [1,2,3,4,5]. The hydraulic conductivity, a key parameter in seepage analysis, is closely related to the accuracy of the analysis results. Minimized error between the simulated and actual values of hydraulic conductivity could improve the reliability of the analysis. In-situ tests have been proven to be helpful in determining the hydraulic conductivity of dams and dam foundations. The hydraulic conductivity determined by in-situ testing agrees with the actual value when the test samples are small. However, this method is time-consuming and costly when there are large quantities of models. Another method to solve this problem is inverse analysis. The inverse research based on monitoring data and numerical simulation results demonstrates economy and efficiency. The essence of the inverse analysis is to determine the hydraulic conductivity by measurements and simulated results. In the inversion analysis, optimization algorithms are widely applied for iterative search over a range of hydraulic conductivity values. The optimal hydraulic conductivity is determined by iteration while minimizing the objective function.

Considering the repetition of the iterative process, optimization algorithms are widely used to improve efficiency and accuracy in the searching process of the inverse problem. For example, on the basis of the Radial Basis Function (RBF) neural network optimized by Particle Swarm Optimization (PSO), Chi et al. [6] constructed an inverse model for the permeability coefficient of a high core rockfill dam; Combining error Back-Propagation Neural Network (BPNN) and Genetic Algorithm (GA), Deng and Lee [7] proposed an inverse analysis method for determining the displacements. This method was successfully applied in the displacement identification of the lock profile of the Three Gorges Project, which led to reasonable results. Zhao et al. [8] developed the differential evolution (DE) algorithm to determine soil parameters in the field of deep excavation, which improved the stability of the backtracking parameters. Simulated annealing [9,10,11,12] and ant colony optimization [13,14,15,16,17] have also been extensively used in the inverse problem. Significantly, much progress has been made in the research field of seepage because of optimization algorithms [18,19,20,21,22,23,24,25]. Tan et al. [26] proposed a biological immune mechanism-based quantum particle swarm optimization (IQPSO) algorithm to solve the inversion problem of seepage parameters. Based on back propagation neural network (BPNN) and genetic algorithm (GA), Zhou et al. [27] developed a new approach for inverse modeling of the transient groundwater flow in dam foundations, which improved the uniqueness and reliability of the inversed results and made tractable the large-scale inverse problems in engineering practices. Zhang et al. [28] proposed an inverse analysis model by using the genetic algorithm (GA) and finite element analysis technology, to solve the calcium leaching problems.

Although optimization algorithms are frequently employed for inverse problems, they suffer from low accuracy, slow convergence, and poor robustness. The Gray Wolf Optimizer (GWO) proposed by Mirjalili [29] has been shown to be efficient and intelligent in engineering optimization. Mirjalili [29] compared the performance of GWO with Particle Swarm Optimization (PSO), Gravity Search Algorithm (GSA), Differential Evolution (DE), Evolutionary Programming (EP), and Evolutionary strategy (ES). The results demonstrate that the GWO can provide very competitive results compared to these well-known meta-heuristics. It has been extensively adopted in various fields due to its simple structure, fewer parameters, and easy coding implementation. However, it tends to converge to locally optimal solutions. In addition, suboptimal values could result from completely randomized initial populations. Therefore, strategies of improvement are proposed as necessary.

Generally, there are three main strategies to improve the GWO, including adjustments of initial populations, convergence factor, and formula of a location update [30]. Pradhan et al. [31] combined the concept of opposition with GWO, initially providing a uniform population for the algorithm. Long et al. [32] introduced the theory of good point set to population initialization, which improved the homogeneity of the population. Based on this, Long et al. [33] considered the dynamics of the iterative process and proposed an equation of the convergence factor based on the number of iterations. This exponential function simulates the iterative process and balances the local and global search to a certain extent. Mittal et al. [34] described the decay process of the parameter by an exponential function, which leads to improved accuracy for GWO. Salgotra et al. [35] applied the spiral property from the whale optimization algorithm (WOA) to the GWO, which solves the premature convergence in the evolutionary algorithm. Mostafa et al. [36] introduced variational operators to update the location of the individual in GWO and improved the algorithm’s performance. Gupta and Deep [37] adopted a random wandering strategy to enhance the accuracy of the algorithm.

A great deal of research has been conducted on the improvement of GWO. However, local and global search, along with homogenization and randomization, cannot be relatively balanced by these improvements. GWO based on these strategies is still limited in terms of efficiency and accuracy. Therefore, there remains potential for improvement across these dimensions.

Based on the evolution of the GWO, three strategies are proposed to ensure the accuracy and efficiency of the algorithm. Initial populations of semi-uniform and semi-random were proposed for rational initialization. A polynomial-based nonlinear convergence factor was applied to maintain a balance between global and local search. Weighting factors based on Euclidean norms and hierarchy were given to dynamically update the wolves’ positions for jumping out of the local optimum at the late stage of the iterative search [38,39]. The Improved Gray Wolf Optimizer (IGWO) has been proven to be effective through numerical experiments. This algorithm was applied for the inversion model of the Kakiwa Dam located in Sichuan Province, China. The fitness function is constructed by measurements of piezometers and measuring weir. The finite element method was used to simulate the seepage process under the assumption of steady flow. The free tetrahedral grid is used to construct the two-dimensional mesh of the dam. The finite element calculation is carried out in COMSOL Multiphysics, while the iterative control and data extraction are implemented in MATLAB. The objective hydraulic conductivity and the corresponding fitness value were determined when the maximum number of iterations was achieved. The seepage field at the dam site was also presented.

2. Improved Gray Wolf Optimizer

2.1. Overview of Gray Wolf Optimizer

The Gray Wolf Optimizer (GWO) is a new group intelligence algorithm considering gray wolves’ hierarchy and group hunting. There is a strict hierarchy in the gray wolf population. The population is classified into four levels of status in accordance with the fitness values of individuals. Wolves in the first level of the population are responsible for making decisions and leading the group in hunting. Wolves in the second level take responsibility for helping to manage the group. Wolves in the third level of the group obey the orders of the first two levels of gray wolves. All the remaining populations are set at the fourth level.

The hunting process tends to be taken as group action in gray wolf populations, which could be summarized by tracking, encircling, and attacking. The rank in the wolf pack changes dynamically with the individual fitness value in the hunting process. The fitness value can be considered the distance between a wolf and its prey. This means that the closer the distance to the target, the higher the level of the wolf.

Let $\alpha$, $\beta$, and $\gamma$ represent the three dominant wolves in rank order. The mathematical model of gray wolf hunting is established as described in Equation (1).

(1)$\{\begin{array}{c}d=\left|C·X_p\left(t\right)-X\left(t\right)\right|\\ X\left(t+1\right)=X\left(t\right)-A·d\end{array}$

where $d$ is the perceived distance between the gray wolf and the prey; $t$ represents the number of iterations. $C$ denotes the coefficient vector; $X_p\left(t\right)$ and $X\left(t\right)$ are the positions of the wolf and the prey, respectively; $X\left(t+1\right)$ means the position of the wolf after the iteration; $A$ stands for the coefficient vector.

The expressions for $A$ and $C$ are shown in Equation (2).

(2)$\{\begin{array}{c}A=a·\left(2m_1-1\right)\\ C=2m_2\end{array}$

where $a$ is the convergence factor; $m_1$ and $m_2$ are both random numbers between 0 and 1.

It is assumed that the first three levels of gray wolves have a better perception of the location of the prey. The populations in the fourth layer decide the direction and distance of the next movement according to the positions of the three dominant wolves. The wolf’s position in the fourth rank is updated according to Equation (3).

(3)$X\left(t+1\right)=\frac{1}{3}{\displaystyle \sum }_{i=1}^3{X}_i$

where $X\left(t+1\right)$ denotes the position of the wolf $\phi$ after update; $X_i\left(i=1,2,3\right)$ represent the position vectors of the wolf $\alpha$, $\beta$, and $\gamma$, respectively.

The position vectors of $\alpha$, $\beta$, and $\gamma$ are expressed by Equation (4).

(4)$\{\begin{array}{c}{X}_1=\left|X_{\alpha }-A_1·d_1\right|\\ X_2=\left|X_{\beta }-A_2·d_2\right|\\ X_3=\left|X_{\gamma }-A_3·d_3\right|\end{array}$

where $X_i\left(i=1,2,3\right)$ stand for the position vectors of the wolf $\alpha$, $\beta$, and $\gamma$, respectively; $X_j\left(j=\alpha ,\beta ,\gamma \right)$ represent the prey position perceived by $\alpha$, $\beta$, and $\gamma$, respectively. $A_i\left(i=1,2,3\right)$ are the coefficient vectors; $d_i\left(i=1,2,3\right)$ mean the distance between the three dominant wolves and the prey.

$d_i\left(i=1,2,3\right)$ could be expressed as Equation (5) follows.

(5)$\{\begin{array}{c}{d}_1=\left|C_1·X_{\alpha }-X\left(t\right)\right|\\ d_2=\left|C_2·X_{\beta }-X\left(t\right)\right|\\ d_3=\left|C_3·X_{\gamma }-X\left(t\right)\right|\end{array}$

where $C_i\left(i=1,2,3\right)$ are the coefficient vectors; $X\left(t\right)$ is the wolf’s position in the t-th iteration.

2.2. Strategies of Improvement

2.2.1. Initial Populations of Semi-Uniform and Semi-Random

The instability of the solution could be increased by completely randomized initial populations, which leads to unstable results. A novel approach for initializing populations was presented to balance uniformity and randomness. The solution range is divided into intervals equidistantly according to the population size. The j-th interval can be expressed by Equation (6).

(6)${\Delta }_j=\left[lb+\frac{ub-lb}{S}\left(j-1\right),lb+\frac{ub-lb}{S}j\right]\left(j=1,2\dots ,S\right)$

where ${\Delta }_j$ is the j-th interval; $lb$ and $ub$ are the upper and lower bounds of the solution set, respectively; $S$ is the population size.

Generate a random initial solution in each interval to ensure randomness. The initial populations are uniformly distributed in the solution space without losing randomness. The expression of the initial solution is given in Equation (7).

(7)$R_j=lb+\frac{ub-lb}{S}j+rand()·\frac{ub-lb}{S}\left(j=1,2\dots ,S\right)$

where $R_j$ is the initial solution of the j-th interval; $lb$ and $ub$ are the upper and lower bounds of the solution set, respectively; $S$ is the population size; $rand()$ is a random real number between 0 and 1.

2.2.2. Polynomial-Based Nonlinear Convergence Factor

The way to search for prey is determined by the coefficient vector $A$. The gray wolf can be in any position between the current individual and the prey at the next moment when $\left|A\right|<1$, indicating that the next position of the wolf will be closer to the location of the prey. This is considered a local search. When $\left|A\right|>1$, the next location of the wolf will be further away from the prey than the current location. Gray wolves tend to search over a wider area, which is considered a global search. The positions of gray wolves change with a rapid convergence speed. In this case, the search step size of the gray wolf becomes smaller, thus achieving a refined search. $A$ varies dynamically with convergence factor $a$. Considering good symmetry and smoothness, a function based on a third-degree polynomial was used to fit the convergence factor $a$. The expression for $a$ is shown in Equation (8).

(8)$a\left(t\right)=r_1{(\frac{t}{T})}^3+r_2{(\frac{t}{T})}^2+r_3\frac{t}{T}+r_4$

where $a$ is the convergence factor; $r_i\left(i=1,2,3,4\right)$ denote the real-valued parameters. $t$ represents the number of iterations. $T$ is the maximum number of iterations.

The constraints are given here, as presented in Equation (9).

(9)$\{\begin{array}{c}a\left(0\right)=2,a\left(T\right)=0,a\left(\frac{T}{2}\right)=1\\ a^{\prime }\left(0\right)<0,a^{\prime }\left(T\right)<0\end{array}$

Thus, $a$ could be indicated by Equation (10).

(10)$$\begin{gathered} a\left(t\right)=\frac{4+2r_3}{T^3}{t}^3-\frac{6+3r_3}{T^2}{t}^2+r_3\frac{t}{T}+2 \\ -3.0\times {10}^{-3}<r_3<0 \end{gathered}$$

Figure 1 shows the nonlinear convergence factor evolution at different values of $r_3$. The maximum number of iterations is set to 500. The three values of $r_3$ are $-2\times {10}^{-3}$, $-1\times {10}^{-3}$, and $-2\times {10}^{-4}$. The convergence factor values for the three curves decrease with the number of iterations. With the increase of iterations, the cut-off point between global and local search is reached when $T_0=250$. Global and local searches could be equally divided and effectively balanced in this condition. In addition, the curve corresponding to $r_3=-2\times {10}^{-3}$ is lower than the other two curves at the early search stage and higher at the late search, indicating its focus on local search and adequate step size. Considering the drawback of converging to the local optimum in GWO, the curve corresponding to $r_3=-2\times {10}^{-3}$ was chosen to ensure the property of jumping out of the local optimum in the study.

2.2.3. Weighting Factors Based on Euclidean Norm and Hierarchy

In the GWO strategy, the position of the gray wolf is updated by the average formula. One limitation of this strategy is that the leadership of the wolf located in the first rank is not considered. Another weakness is that the weighting factors are kept constant as the iteration proceeds. Weighting factors based on Euclidean norms and hierarchy are proposed to overcome this problem, as shown in Equation (11).

(11)$\{\begin{array}{c}{\rho }_1=\frac{\Vert X_1\Vert }{\Vert X_1\Vert +\Vert X_2\Vert +\Vert X_3\Vert }\\ {\rho }_2=\frac{\Vert X_2\Vert }{\Vert X_1\Vert +\Vert X_2\Vert +\Vert X_3\Vert }\\ {\rho }_3=\frac{\Vert X_3\Vert }{\Vert X_1\Vert +\Vert X_2\Vert +\Vert X_3\Vert }\end{array}$

Here ${\rho }_i \left(i=1,2, 3\right)$ are the weighting factors of the first-three level wolves, respectively; $X_i \left(i=1, 2, 3\right)$ are the Euclidean norms of the first three levels, respectively.

The weights of the wolf $\alpha$, $\beta$, and $\gamma$ are multiplied by 0.6, 0.3, and 0.1 to reinforce the wolf pack hierarchy. The formula for updating the location of the gray wolf can be improved, as Equation (12) states.

(12)$X\left(t+1\right)=\frac{6{\rho }_1·X_1+3{\rho }_2·X_2+{\rho }_3·X_3}{10}$

where ${\rho }_i \left(i=1,2, 3\right)$ are the weighting factors of the first-three level wolves, respectively; $X_i \left(i=1, 2 , 3\right)$ represent the position vectors of the wolf $\alpha$, $\beta$, and $\gamma$ respectively.

2.3. Numerical Experiment of Algorithm Performance

A numerical experiment was performed to verify the effectiveness of IGWO. Six typical functions were selected for simulation in the experiments, including Sphere, Rosenbrock, Quartic, Rastrigin, Ackley, and Griewank. Table 1 shows the mathematical expressions, dimensions, and search ranges of these typical functions. Sphere, Rosenbrock, and Quartic are single-peak functions. Especially Quartic is a multidimensional flat bottom function with random disturbances. The single-peak functions are mainly applied to determine the accuracy of IGWO. Rastrigin, Ackley, and Griewank are multi-peaked functions that tend to cause the algorithm to converge to a locally optimal solution. The performance to jump out of the local optimum could be tested reasonably for IGWO. In addition, the experimental results of IGWO, SGWO [33], and GWO are compared in the simulation.

To ensure a fairness, the population size is 30, and the maximum number of iterations is 500 for IGWO, SGWO, and GWO. The three algorithms were performed 30 times independently for each function, and the average values were taken as the simulation results. The results of the numerical experiments are given in Table 2. The optimal values simulated by IGWO are closer to the theoretical optimal solution than the results of the other two algorithms. Among them, the simulation results of Sphere, Rastrigin, and Griewank are equal to the theoretical values, revealing the high accuracy of IGWO.

Figure 2 presents the convergence curves of the six test functions with IGWO, SGWO, and GWO. The convergence curves of all three algorithms continue to decrease as the number of iterations increases. In particular, the convergence curve of IGWO decreases significantly faster than the corresponding curves of the other two algorithms, indicating the progress of IGWO in terms of running time. Furthermore, for the multi-peaked test function, the convergence curve of IGWO continues to decrease while the search of the other two algorithms converges to a local optimum solution. This result demonstrates the ability to jump out of the local search and converge to the global optimum value in IGWO.

3. Inverse Model of Seepage Parameters

3.1. The Objective Function

The aim of the inverse model is to determine the hydraulic conductivity for each partition while minimizing the value of the objective function. The objective function was constructed by hydraulic head and leakage in this paper, improving the reliability of the simulation results. Suppose that the hydraulic conductivity of each medium is isotropic. Denote by $\mathit{K}=\left[k_1, k_2,\dots , k_n\right]$ the combination of hydraulic conductivity, in which $k_i$ represents the hydraulic conductivity of the ith media. The number of piezometers and measuring weirs are indicated by $m$ and $n$, respectively. $\mathit{H}=\left[H_1, H_2,\dots , H_M\right]$ is expressed as a sequence of hydraulic head measurements. Similarly, $\mathit{Q}=\left[Q_1, Q_2,\dots , Q_N\right]$ is a series of leakage volume measurements. $H_i\left(\mathit{K}\right)$ and $Q_j\left(\mathit{K}\right)$ are the simulated hydraulic head and leakage volume by finite element method. The mathematical model for the inverse problem is established, as shown in Equation (13).

(13)$$\begin{gathered} minf={\left({\displaystyle \sum }_{i=1}^M\frac{\Vert H_i\left(\mathit{K}\right)-H_i{\Vert }_2^2}{\Vert H_i\Vert {}_2^2}\right)}^{\frac{1}{2}}+w{\left({\displaystyle \sum }_{i=1}^N\frac{\Vert Q_j\left(\mathit{K}\right)-Q_j\Vert {}_2^2}{\Vert Q_j\Vert {}_2^2}\right)}^{\frac{1}{2}} \\ \mathrm{s}.\mathrm{t}. {\mathit{K}}_{min}\le \mathit{K}\le {\mathit{K}}_{max} \end{gathered}$$

Here ${\mathit{K}}_{min}$ and ${\mathit{K}}_{max}$ are the lower and upper bounds of hydraulic conductivity values, respectively. The range of hydraulic conductivity can be roughly determined by geological data and engineering experience. $w$ is a weight factor for balancing the hydraulic head and the leakage volume. In this paper, the simulated leakage value is estimated by the flow rate and area of the overwater cross-section. It is suggested that the value of the weighting factor is set small considering an error between the simulated value and the measurement at the shoulder part of the dam. Zhou [23] compared the relative errors of hydraulic head and leakage volume at different weights. The results show that the simulated values are in good agreement with the measurements, and the minimum value of relative error is reached at the condition of $w=0.02$. The finding was applied in this paper.

The objective combination of hydraulic conductivity was obtained by the searching process. The search process was accelerated by IGWO. The objective combination of hydraulic conductivity was applied to simulate the seepage field of the dam during operation.

3.2. Procedure of the Inversion Model

The procedure of the inversion model could be summarized by specific steps. The steps are as follows.

Step 1: Set initial parameters of IGWO. The number of search agents $S$, the maximum number of iterations $T$, and the bounds of hydraulic conductivity values, ${\mathit{K}}_{min}$ and ${\mathit{K}}_{max}$, are determined initially.

Step 2: Initialize the population. Equation (1) is applied for initialization, ensuring that the initial populations are uniformly distributed in the solution space.

Step 3: Calculate the fitness of individual gray wolves. The three gray wolves with the top fitness values are selected as $\alpha$, $\beta$, and $\gamma$.

Step 4: Update the position. Determine the distance of the gray wolf from the three dominant wolves, respectively. Orient the location and calculate the weighting factors of the three wolves. The position of the gray wolf at the fourth level is updated according to Equation (12).

Step 5: Iterative Judgment. Determines whether the maximum number of iterations has been achieved. If not, skip to the third step. Otherwise, end the iterative procedure. Output the objective hydraulic conductivity and the corresponding fitness value.

Step 6: Positive verification. The combination of the target hydraulic conductivity is substituted into the finite element model for positive verification. Compare the calculated and monitored values of hydraulic head and leakage volume and evaluate the reasonableness of the simulation.

Figure 3 presents the flow chart of the model.

4. A Casebook Study

4.1. Project Overview

Located on the Muli River in the Sichuan Province of China, the Kakiwa Hydropower Station is a project focused on power generation and ecological preservation. The location of the Kakiwa Dam is indicated in Figure 4.

A concrete panel rockfill dam is selected as the barrage in the pivot project, with a maximum height of 171 m. The crest width of the dam is 11 m. The dam mainly consists of a concrete face slab, blanket area, cushion area, transition area, rockfill area, drainage area, ballast area, and grout curtain. The normal storage level is 2850.00 m, the calibration flood level is 2852.20 m, and the dead water level is 2800.00 m. Figure 5 shows the maximum cross-section of the dam body.

A total of 37 piezometers were installed for seepage monitoring of the dam body and dam foundation. Among them, the piezometer $P_{DB-13}$ is located downstream of the curtain. $P_{DB-24}-P_{DB-27}$ are installed near the original ground line. Figure 5 presents the locations and elevations of these piezometers. The water measuring weir is installed at the downstream cofferdam axis of the dam.

The hydraulic conductivities of the five media, including the grout curtain, the top cover layer, the second cover layer, the moderately weathered zone, the slightly weathered zone, and the fresh bedrock zone, are limited to reasonable ranges and needed to be optimized. The parameter ranges are given in the results of the simulation. In addition, the hydraulic conductivities of other media are indicated in Table 3.

4.2. Analysis of Monitoring Data

Figure 6 shows the monitoring data of the piezometers around the grout curtain and the water measuring weirs during the operation period. The hydraulic head on the downstream side of the curtain is relatively consistent with the upstream water level. The hydraulic head measured by the piezometer lags behind the upstream head, which is called the hysteresis effect. The measured hydraulic head rises less than the upstream water level. In addition, the value of the piezometer $P_{DB-24}$ is approximately 2776 m, and the difference among the values of piezometers $P_{DB-24}$ to $P_{DB-27}$ is not significant, indicating the efficiency of the impermeable curtain.

Similarly, the seepage volume is consistent with the upstream water level. The value of seepage volume at the dam body is relatively small, with a stable value of 0.41 m³/s in the operation period.

It is assumed that the seepage field is stable for simplicity. A period with a slight variation of upstream and downstream reservoir levels and long duration was selected for the inverse model, which could minimize the seepage lag effect to a certain extent. As seen in Figure 6, the upstream pool level is relatively stable from 10 July 2017, to 27 October 2017, with values between 2849.55 m and 2850.06 m. Therefore, this period was chosen for the inversion model.

4.3. Computation Model

IGWO was used for iterative search. The maximum number of iterations is 200, and the initial population size is 30. The simulation is performed 20 times independently by IGWO, and the average values of simulated hydraulic conductivity are taken as the target results.

The multi-physics field simulation software COMSOL Multiphysics was applied in this simulation. A two-dimensional finite element mesh was created for inverse modeling each medium’s hydraulic conductivity. The seepage properties of the dam body and foundation were analyzed. The free tetrahedral grid is applied to build the mesh, composed of 11,635 domain elements and 14,400 vertices. The maximum and minimum widths of the cells are 20 m and 5 m, respectively. A right-hand cartesian coordinate system is constructed with the x-axis pointing to the downstream reservoir and the z-axis pointing to the sky vertically. The upstream and downstream are both extended by 2.5 times the dam’s height on the x-axis. The depth of the foundation is taken as 400 m. In addition, the upstream and downstream water levels are, respectively, set at 2850 m and 2702 m. The two-dimensional finite element mesh is presented in Figure 7.

4.4. Results of the Simulation

4.4.1. Hydraulic Conductivity

The hydraulic conductivity of each medium was determined by IGWO. The hydraulic conductivity of each medium at the dam site is given in Table 4. All the results are within the corresponding search range.

4.4.2. Hydraulic Head

The results of hydraulic conductivity were substituted into the finite element model for positive analysis to verify the reasonableness of this simulation. The simulated values of the hydraulic head and leakage volume at the monitoring points were obtained and compared with the corresponding measurements. Absolute and relative errors of hydraulic head and leakage were calculated. A contour of the hydraulic head in the dam site area was also predicted.

The mathematical expression for the relative error of the hydraulic head is given in Equation (14):

(14)${\delta }_H=\frac{\left|H_i-H\right|}{\Delta H}\times 100\%$

where ${\delta }_H$ is the relative error of the hydraulic head; $H_i$ and $H$ represent the simulated and measured values of the hydraulic head, respectively. $\mathsf{\Delta }H$ denotes the difference in water level between the upstream and downstream sides, taken as 148 m.

Table 5 compares the calculated and measured hydraulic head values at the monitoring points. The calculated values at the monitoring points are relatively close to the measurements. Among them, the simulated values at the piezometers $P_{DB-24}-P_{DB-27}$ show a very high consistency with the corresponding measurements. The maximum value of the absolute error −1.41 m, and the maximum value of the relative error is 0.95%. Meanwhile, the value of the hydraulic head decreases with the increase of seepage distance.

4.4.3. Leakage Volume of the Dam Foundation

The seepage volume is estimated by the flow velocity and the overflow surface. The formula for calculating leakage volume is given in Equation (15).

(15)$Q_i={\displaystyle \iint }vBdxdz$

where $Q_i$ is the simulated seepage volume; $v$ means the flow rate at the vertical spillway surface; $B$ represents the length of the dam taken as 355 m.

The relative error for the leakage volume is determined by Equation (16).

(16)${\delta }_Q=\frac{\left|Q_i-Q\right|}{Q}\times 100\%$

where ${\delta }_Q$ is the relative error of the dam body leakage; $Q_i$ and $Q$ are the simulated and measured values of the dam body leakage, respectively.

Table 6 shows the comparison of the calculated and measured values of the dam body leakage during the stable upstream water level. The simulated leakage values are in good agreement with the actual measurements, showing the accuracy of the IGWO strategy and the reliability of the simulation.

The leakage measurements are averaged over the simulation period. The relative error of the calculated seepage volume is 5.00%, demonstrating the positive performance of IGWO in the simulation.

Figure 8 presents the contour map of the hydraulic head at the dam site. The distribution of contours is in accordance with seepage characteristics. The results show the accuracy of the simulation and the reasonableness of the calculated hydraulic conductivity.

5. Discussion

The hydraulic conductivities of the dam and dam foundation change with operation and loading, affecting the effectiveness and safety of the hydraulic project. An Improved Gray Wolf Optimizer for solving the problem that hydraulic conductivity is not easily determined, is introduced in this paper.

The objective hydraulic conductivity was determined based on IGWO. The obtained hydraulic conductivities of the dam and dam foundation were applied in the finite element positive analysis. The results show the effectiveness of IGWO in determining hydraulic conductivity. It is suggested that IGWO could be used to obtain reasonable simulation results in similar inverse problems.

The hydraulic conductivity of each medium in the dam body and dam foundation simulated in this paper was obtained with a stable upstream water level. In fact, the upstream water level of a reservoir varies continuously with the operating conditions and purpose. The hydraulic conductivity will also change with the monitored values. It is worth noting that IGWO is still applicable under this condition.

6. Conclusions

Aiming at balancing the local and global search of GWO, along with uniformity and stochasticity, IGWO was proposed in this paper. The improvement of IGWO in accuracy and running time was indicated in numerical experiments. IGWO was used in the inverse modeling of the hydraulic conductivity of the Kakiwa dam. The hydraulic head and leakage were used to set up the objective function. Errors in the hydraulic head and leakage were calculated. The steady seepage field was analyzed in the application case. The main conclusions of this paper are given as follows.

(1) The performance of IGWO is improved due to the three strategies. A novel approach for initialization contributes to populations of semi-uniform and semi-random. The polynomial-based nonlinear convergence factor is selected to keep the equilibrium of the local and global search. Weighting factors based on Euclidean norms and hierarchy helps to update the position of the wolf dynamically.

(2) A numerical experiment was conducted to demonstrate the performance of IGWO. The optimal values obtained by IGWO are closer to the theoretical solution than the results of the other two algorithms. The results of the experiment demonstrate the feasibility and efficiency of IGWO.

(3) The target combination of hydraulic conductivity was obtained by IGWO. The values of hydraulic conductivities were substituted into the finite element model. The values of the hydraulic head and leakage quantity at the corresponding measurement points were obtained. The maximum values of absolute and relative errors of the hydraulic head were—1.41 m and 0.95%, respectively. The absolute and relative errors of the seepage volume were—0.02 m³/s and 5.00%, respectively. The results of the application case show that the inversion model and the algorithm are reliable and efficient.

Author Contributions

Conceptualization, Y.S. and L.X.; methodology, Z.S.; software, Y.S.; validation, J.D. and L.J.; formal analysis, L.X.; investigation, J.D.; resources, Z.S.; data curation, L.J. and Q.L.; writing—original draft preparation, Y.S.; writing—review and editing, Y.S.; visualization, J.D. and Q.L.; supervision, Z.S.; project administration, Z.S.; funding acquisition, L.X. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Nonlinear convergence factor based on third-degree polynomials.

Figure 2. Convergence curves of the test function with the three optimization algorithms: (a) Sphere; (b) Rosenbrock; (c) Quartic; (d) Rastrigin; (e) Ackley; (f) Griewank.

Figure 2. Convergence curves of the test function with the three optimization algorithms: (a) Sphere; (b) Rosenbrock; (c) Quartic; (d) Rastrigin; (e) Ackley; (f) Griewank.

Figure 3. Flow chart of the inverse model based on IGWO.

Figure 4. The location of the Kakiwa Dam.

Figure 5. The max cross-section of the dam body.

Figure 6. Measurements of hydraulic head and seepage discharge during storage period: (a) hydraulic head; (b) seepage volume.

Figure 7. 2D FE mesh of Kakiwa dam.

Figure 8. Distribution of seepage contours at the maximum section of the dam.

Table 1

Test functions in the numerical experiment.

Test Function	Mathematical Expression	Dimension	Search Range
Sphere	${\displaystyle {\displaystyle \sum }_{i=1}^D{x}_i^2}$	30	[–100, 100]
Rosenbrock	${\displaystyle {\displaystyle \sum }_{i=1}^{D-1}\left[100{\left(x_{i+1}-x_i^2\right)}^2+{\left(x_i-1\right)}^2\right]}$	30	[−30, 30]
Quartic	${\displaystyle {\displaystyle \sum }_{i=1}^{D}i{x}_i^4+random\left[0, 1\right]}$	30	[−1.28, 1.28]
Rastrigin	${\displaystyle {\displaystyle \sum }_{i=1}^{D-1}\left[x_i^2-10cos\left(2\pi x_i\right)+10\right]}$	30	[−5.12, 5.12]
Ackley	${\displaystyle -20exp\left(-0.2\sqrt{\frac{1}{30}{\displaystyle \sum }_{i=1}^D{x}_i^2}\right)-exp\left(\frac{1}{30}{\displaystyle \sum }_{i=1}^{D}cos2\pi x_i\right)+20+e}$	30	[−32, 32]
Griewank	${\displaystyle \frac{1}{4000}{\displaystyle \sum }_{i=1}^D{x}_i^2-{\displaystyle \prod }_{i=1}^{D}cos\left(\frac{x_i}{\sqrt{i}}\right)+1}$	30	[−600, 600]

Table 2

Results of the numerical experiments.

Test Function	Optimization Algorithm	Simulated Optimum Value	Theoretical Optimum Value
Sphere	IGWO	0	0
	SGWO	6.45 × 10−33
	GWO	1.34 × 10−26
Rosenbrock	IGWO	2.89 × 10	0
	SGWO	2.70 × 10
	GWO	2.72 × 10
Quartic	IGWO	2.04 × 10−4	0
	SGWO	2.83 × 10−4
	GWO	1.40 × 10−3
Rastrigin	IGWO	0	0
	SGWO	5.68 × 10−14
	GWO	1.71 × 10−12
Ackley	IGWO	4.44 × 10−15	0
	SGWO	1.51 × 10−14
	GWO	1.11 × 10−13
Griewank	IGWO	0	0
	SGWO	1.16 × 10−2
	GWO	2.84 × 10−2

Table 3

Hydraulic conductivity of stationary medium.

Material	Hydraulic Conductivity(m/s)
Upstream face slab	1.00 × 10−6
Bedding material	9.90 × 10−4
Transition material	7.50 × 10−2
Main rock-fill zone	8.70 × 10−1
Secondary rock-fill zone	9.80 × 10−1
Upstream blanket	1.00 × 10−5

Table 4

Hydraulic conductivity of each medium determined by IGWO.

Material	Hydraulic Conductivity(m/s)	Search Range(m/s)
Grout curtain	3.00 × 10−5	1.00 × 10−6–1.00 × 10−4
The top cover layer	5.33 × 10−2	1.00 × 10−3–1.00 × 10−1
The second cover layer	2.67 × 10−3	1.00 × 10−4–1.00 × 10−2
Moderately-weathered zone	5.50 × 10−4	1.00 × 10−5–1.00 × 10−3
Slightly weathered zone	1.10 × 10−4	1.00 × 10−5–1.00 × 10−3
Fresh bedrock zone	3.10 × 10−5	1.00 × 10−6–1.00 × 10−4

Table 5

Comparison between the calculated and measured hydraulic head values.

MonitoringPoints	MeasuredHydraulic Head(m)	Simulated Hydraulic Head(m)	Absolute Error(m)	Relative Error(%)
$P_{DB-13}$	2736.17	2735.74	−0.43	0.29
$P_{DB-24}$	2718.67	2717.51	−1.16	0.78
$P_{DB-25}$	2718.28	2717.16	−1.12	0.76
$P_{DB-26}$	2718.21	2716.80	−1.41	0.95
$P_{DB-27}$	2716.61	2716.24	−0.37	0.25

Table 6

Comparison between the simulated and measured values of the dam body leakage.

Leakage	Measured Values (m3/s)	Simulated Values(m3/s)	Absolute Error(m3/s)	Relative Error(%)
Dam body	0.40	0.38	−0.02	5.00

References

1. Liu, M.; Zhong, J.Q.; Rao, X.B.; Liu, J.; Wang, W.J. Seepage monitoring analysis and safety assessment for Yatang dam. J. Yangtze River Sci. Res. Inst.; 2007; 24, pp. 30-33. (In Chinese)

2. Chu, X.W.; Xu, M.; Wang, Z.M. Stability analysis of seepage on the accumulation dam of a phosphogypsum tailings. J. Eng. Geol.; 2016; 24, pp. 661-667. (In Chinese)

3. Xu, J.B.; Wei, W.; Bao, H.; Zhang, K.; Lan, H.; Yan, C.; Sun, W. Failure models of a loess stacked dam: A case study in the ansai area (China). Bull. Eng. Geol. Environ.; 2020; 79, pp. 1009-1021. [DOI: https://dx.doi.org/10.1007/s10064-019-01605-z]

4. Jiang, Z.X.; He, J.P. Detection model for seepage behavior of earth dams based on data mining. Math. Probl. Eng.; 2018; 2018, 8191802. [DOI: https://dx.doi.org/10.1155/2018/8191802]

5. Fan, X.F.; Wu, Z.Y.; Liu, L.; Wen, Y.; Yu, S.; Zepeng, Z.; Li, Z. Analysis of sluice foundation seepage using monitoring data and numerical simulation. Adv. Civ. Eng.; 2019; 2019, 2850916.

6. Chi, S.; Ni, S.; Liu, Z. Back analysis of the permeability coefficient of a high core rockfill dam based on a RBF neural network optimized using the PSO algorithm. Math. Probl. Eng.; 2015; 2015, pp. 1-15. [DOI: https://dx.doi.org/10.1155/2015/124042]

7. Deng, J.H.; Lee, C.G. Displacement back analysis for a steep slope at the Three Gorges Project site. Int. J. Rock Mech. Min.; 2001; 38, pp. 259-268. [DOI: https://dx.doi.org/10.1016/S1365-1609(00)00077-0]

8. Zhao, B.D.; Zhang, L.L.; Jeng, D.S.; Wang, J.H.; Chen, J.J. Inverse analysis of deep excavation using differential evolution algorithm. Int. J. Numer. Anal. Met.; 2015; 39, pp. 115-134. [DOI: https://dx.doi.org/10.1002/nag.2287]

9. Ortiz-Aleman, C.; Martin, R. Inversion of electrical capacitance tomography data by simulated annealing: Application to real two-phase gas-oil flow imaging. Flow Meas. Instrum.; 2005; 16, pp. 157-162. [DOI: https://dx.doi.org/10.1016/j.flowmeasinst.2005.02.014]

10. Tabatabaeenejad, A.; Moghaddam, M. Inversion of subsurface properties of layered dielectric structures with random slightly rough interfaces using the method of simulated annealing. IEEE Trans. Geosci. Remote Sens.; 2009; 47, pp. 2035-2046. [DOI: https://dx.doi.org/10.1109/TGRS.2008.2011982]

11. Song, G.R.; Hong, G.F.; Lu, Y.; Xu, Y.; Qin, D.; Wu, B.; He, C. Inversion of elastic constants of anisotropic (100) silicon based on surface wave velocity by acoustic microscopy using Particle Swarm-Based-Simulated annealing optimization. J. Nondestruct. Eval.; 2015; 34, 43. [DOI: https://dx.doi.org/10.1007/s10921-015-0316-6]

12. Netto, A.; Dunbar, J. 3-D constrained inversion of gravimetric data to map the tectonic terranes of southeastern Laurentia using simulated annealing. Earth Planet. Sc. Lett.; 2019; 513, pp. 12-19. [DOI: https://dx.doi.org/10.1016/j.epsl.2019.02.018]

13. Shi-chang, L.; Qing-sheng, Z.; Zhe, Y.; Hao-lan, Y. Ant colony optimization for nonlinear AVO inversion of network traffic allocation optimization. Expert Syst. Appl.; 2010; 37, pp. 8343-8347. [DOI: https://dx.doi.org/10.1016/j.eswa.2010.05.049]

14. Lin, P.; Liu, X.L.; Chen, H.X.; Kim, J. Ant colony optimization analysis on overall stability of high arch dam basis of field monitoring. Sci. World J.; 2014; 2014, 483243. [DOI: https://dx.doi.org/10.1155/2014/483243] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25025089]

15. Gao, W. Displacement back analysis for underground engineering based on immunized continuous ant colony optimization. Eng. Optimiz.; 2016; 48, pp. 868-882. [DOI: https://dx.doi.org/10.1080/0305215X.2015.1061814]

16. Ghorbani, E.; Moosavi, M.; Hossaini, M.F.; Assary, M.; Golabchi, Y. Determination of initial stress state and rock mass deformation modulus at lavarak HEPP by back analysis using ant colony optimization and multivariable regression analysis. Bull. Eng. Geol. Environ.; 2021; 80, pp. 429-442. [DOI: https://dx.doi.org/10.1007/s10064-020-01936-2]

17. Li, Z.X.; Rao, S.W. Frequency domain soil parameters inversion of horizontally multilayered earth model with considering high-frequency field. IET Gener. Transm. Dis.; 2018; 12, pp. 5690-5699. [DOI: https://dx.doi.org/10.1049/iet-gtd.2018.5599]

18. He, J.; Chen, S.H.; Shahrour, I. A revised solution of equivalent permeability tensor for discontinuous fractures. J. Hydrodyn.; 2012; 24, pp. 711-717. [DOI: https://dx.doi.org/10.1016/S1001-6058(11)60295-3]

19. Virbulis, J.; Bethers, U.; Saks, T.; Sennikovs, J.; Timuhins, A. Hydrogeological model of the baltic artesian basin. Hydrogeol. J.; 2013; 21, pp. 845-862. [DOI: https://dx.doi.org/10.1007/s10040-013-0970-7]

20. Yu, H.; Wang, X.; Ren, B.; Zeng, T.; Lv, M.; Wang, C. An efficient bayesian inversion method for seepage parameters using a data-driven error model and an ensemble of surrogates considering the interactions between prediction performance indicators. J. Hydrol.; 2022; 604, 127235. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.127235]

21. Xu, L.Q.; Li, L. Inversion Analysis of seepage parameters based on improved shuffled frog leaping algorithm. Math. Probl. Eng.; 2021; 2021, 6536294. [DOI: https://dx.doi.org/10.1155/2021/6536294]

22. Shi, Z.; Wu, Z.; Gu, C.; Chen, B.; Zhang, H.; Yin, W.; Wu, B. Calculation methods for the permeability coefficient of concrete face rockfill dam with cracks. Adv. Civ. Eng.; 2019; 2019, 6571092. [DOI: https://dx.doi.org/10.1155/2019/6571092]

23. Zhang, H.Y.; Song, Z.Y.; Peng, P.; Sun, Y.; Ding, Z.; Zhang, X. Research on seepage field of concrete dam foundation based on artificial neural network. Alex. Eng. J.; 2021; 60, pp. 1-14. [DOI: https://dx.doi.org/10.1016/j.aej.2020.03.041]

24. Kumar, R.; Pal, S.K.; Gupta, P.K. Water seepage mapping in an underground coal-mine barrier using self-potential and electrical resistivity tomography. Mine Water Environ.; 2021; 40, pp. 622-638. [DOI: https://dx.doi.org/10.1007/s10230-021-00788-w]

25. Shu, Y.; Shen, Z.; Xu, L.; Zhang, K.; Yang, C. Inversion analysis of impervious curtain permeability coefficient using calcium leaching model, extreme learning machine, and optimization algorithms. Appl. Sci.; 2022; 12, 3272. [DOI: https://dx.doi.org/10.3390/app12073272]

26. Tan, J.; Xu, L.; Zhang, K.; Yang, C. A biological immune mechanism-based quantum PSO algorithm and its application in back analysis for seepage parameters. Math. Probl. Eng.; 2020; 2020, 2191079. [DOI: https://dx.doi.org/10.1155/2020/2191079]

27. Zhou, C.B.; Liu, W.; Chen, Y.F.; Hu, R.; Wei, K. Inverse modeling of leakage through a rockfill dam foundation during its construction stage using transient flow model, neural network and genetic algorithm. Eng. Geol.; 2015; 187, pp. 183-195. [DOI: https://dx.doi.org/10.1016/j.enggeo.2015.01.008]

28. Zhang, K.; Shen, Z.; Xu, L.; Shu, Y.; Yang, C. Inverse modeling of grout curtain hydraulic conductivity evolution considering the calcium leaching effect. Mathematics; 2022; 10, 381. [DOI: https://dx.doi.org/10.3390/math10030381]

29. Mirjalili, S.; Mirjalili, S.M.; Lewis, A. Grey wolf optimizer. Adv. Eng. Softw.; 2014; 69, pp. 46-61. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2013.12.007]

30. Zhang, X.F.; Wang, X.Y. Comprehensive review of grey wolf optimization algorithm. Comput. Sci.; 2019; 46, pp. 30-38. (In Chinese)

31. Pradhan, M.; Roy, P.K.; Pal, T. Oppositional based grey wolf optimization algorithm for economic dispatch problem of power system. Ain Shams Eng. J.; 2018; 9, pp. 2015-2025. [DOI: https://dx.doi.org/10.1016/j.asej.2016.08.023]

32. Long, W.; Zhao, D.Q.; Xu, S.J. Improved grey wolf optimization algorithm for constrained optimization problem. J. Comput. Appl.; 2015; 35, pp. 2590-2595. (In Chinese)

33. Long, W.; Wu, T.B. Improved grey wolf optimization algorithm coordinating the ability of exploration and exploitation. Control. Decis.; 2017; 32, pp. 1749-1757. (In Chinese)

34. Yin, L.; Sun, Z. Distributed multi-objective grey wolf optimizer for distributed multi-objective economic dispatch of multi-area interconnected power systems. Appl. Soft Comput. J.; 2022; 117, 108345. [DOI: https://dx.doi.org/10.1016/j.asoc.2021.108345]

35. Salgotra, R.; Singh, U.; Singh, S.; Mittal, N. A hybridized multi-algorithm strategy for engineering optimization problems. Knowl.-Based Syst.; 2021; 217, 106790. [DOI: https://dx.doi.org/10.1016/j.knosys.2021.106790]

36. Mostafa, A.; Fouad, A.; Houseni, M.; Allam, N.; Hassanien, A.E.; Hefny, H.; Aslanishvili, I. A hybrid grey wolf based segmentation with statistical image for CT liver images. International Conference on Advanced Intelligent Systems and Informatics; Springer: Cham, Switzerland, 2016; pp. 846-855.

37. Gupta, S.; Deep, K. A novel random walk Grey Wolf Optimizer. Swarm Evol. Comput.; 2019; 44, pp. 101-112. [DOI: https://dx.doi.org/10.1016/j.swevo.2018.01.001]

38. Malik, M.; Mohideen, E.R.; Ali, L. Weighted distance Grey wolf optimizer for global optimization problems. Proceedings of the 2015 IEEE International Conference on Computational Intelligence and Computing Research; Tamilnadu, India, 10–12 December 2015.

39. Panda, M.; Das, B. Grey Wolf Optimizer and Its Applications: A Survey. Proceedings of the 3rd International Conference on Microelectronics, Computing and Communication Systems; Ranchi, India, 12–13 May 2018; Nath, V.; Mandal, J.K. Springer: Singapore, 2019; pp. 179-194.