Inversion of the Permeability Coefficient of a High Core Wall Dam Based on a BP Neural Network and the Marine Predator Algorithm

Abstract

The parameters’ inversion of saturated–unsaturated is important in ensuring the safety of earth dams; many scholars have conducted some research regarding the inversion of hydraulic conductivity based on seepage pressure monitoring data. The van Genuchten model is widely used in saturated–unsaturated seepage analysis, which considers the permeability connected to the water content of the soil and the soil’s shape parameters. A BP neural artificial network is a mature prediction technique based on enough data, and the marine predator algorithm is a new nature-inspired metaheuristic inspired by the movement of animals in the ocean. The BP neural artificial network and marine predator algorithm are applied in the permeability coefficient inversion of a core-rock dam in China; the results show that in the normal operation status, the BP network shows better accuracy, and the average of the absolute error and variance of the absolute error are both minimum values, which are 2.21 m and 1.43 m, respectively. While the water storage speed changes, the marine predator algorithm shows better accuracy; the objective function is calculated to be 0.253. So, the marine predator algorithm is able to accurately reverse the desired results in some situations. According to the actual condition, employing suitable methods for the inverse permeability coefficient of a dam can effectively ensure the safe operation of dams.

1. Introduction

Embankment dams are the most widely used dam style in the world. A core-rockfill dam is a type of embankment dam and is widely used, and core rockfill dams are increasingly being built and put into use. Core-rockfill dams such as Changhe dam (240 m), Nuozhadu dam (261.5 m), and Lianghekou dam (293 m) have been built. Gushui dam (305 m) and Shuangjiangkou dam (314 m) are under construction [1]. Penetration damage is one of the most common forms of damage to embankment dams due to defects in the foundation and body of the dam itself or the effect of dam-penetrating structures, 25% of embankment dam failures [2] are caused by dam failure worldwide, and about 30% of embankment dam wrecks in China are due to penetration damage [3].

In these wrecks, infiltration damage often occurs in the early stage of water storage [4]. When the reservoir first fills with water, the core wall is in a non-saturated state, there is a rapid rise of the reservoir water level, the pore water pressure in the rock pile area increases, the heart wall inside and outside undergo a sudden change in water pressure, and it is prone to infiltration damage. Teton Dam in the United States formed a gradual internal source and developed pipe surges at the contact area of the dam material in the core wall during the initial impoundment of water, leading to breach damage of the dam [5]. An abnormal permeability coefficient of dam-building materials leads to excessive water seepage, which may lead to dam failure. This tragedy might not have happened if the permeability coefficient could have been calculated from monitoring data in time.

The coefficient permeability of unsaturated soils is determined by the water of the soil, which is a function of saturated permeability coefficients and the soil–water characteristic curve (SWCC) [6]. Mualem raised an equation of the relative permeability coefficient $k_r$ from the SWCC [7]. Van Genuchten [8] combined the SWCC formula derived by himself with the Mulam model, raising the specific van-Genuchten model, which is widely used in the calculation of saturated–unsaturated seepage. The Brooks–Corey model [9] is also widely used in the calculation of saturated and unsaturated seepage. Xu [10] used a piece of FEM software named Geo-slope to perform a saturated–unsaturated numerical simulation of a dam during the descent of the reservoir water level; due to the existence of an unsaturated area, the percolation line in the dam is convex and the pore water pressure dissipation becomes slower. Wang [11] simulated slope seepage with variations in the reservoir water level and rainfall intensity based on saturated–unsaturated theory, offering some guidance for slope ability. Shu [12] solved the complicated saturated–unsaturated seepage problem caused by rain infiltration by using the finite element method. Shen [13] built a 3-D FEM calculation model, carrying out a physical model experiment on composite geomembrane defect leakage based on saturated–unsaturated theory.

Though the theory of saturated–unsaturated soil seepage has been developed for a long time, there is still a big problem in calculating the pressure in the dam and slope. In order to accurately determine the pore pressure in the soil, it is necessary to know the seepage calculation parameters of the soil. In saturated–unsaturated seepage analysis, laboratory testing and engineering analogies are the methods for determining the analysis parameters [14,15,16], but these cannot accurately avoid the impact of the size effect, engineering quality, and actual operating environment [17,18,19,20]. Many large-scale projects cannot accurately obtain permeability parameters through experiments, but they can be obtained by inversion through dam monitoring data, while inversion analysis has become an effective method for determining seepage parameters [21].

Many scholars have studied the effect of the water storage rate on the seepage of dams based on saturated–unsaturated theory. Due to the effects of reservoir level fluctuations, there exist unstable seepage flows in the dam body [22]. Yu [23] used the ADALINE network to simulate rainfall infiltration effects and applied it to the processing of dam data. Wang [11] established a diurnal variation monitoring model of dam seepage elements, where a lag effect function based on unstable seepage theory was introduced. Hu [24] used the finite element method to analyze the rainfall and water level changes’ influence on Hualianshu landslide seepage by using the saturated–unsaturated seepage theory. Liu [25] calculated the positional changes of the water table line within the dam body and the corresponding seepage field parameters during the rise in the reservoir water level, discovering the effect of rising levels on upper loose accumulations, which was fairly obvious. Xu [26] simulated the dissipation process of pore water pressure during the rapid decline in the reservoir level and drew relevant security conclusions. Aynaz [27] discovered the relationship between slope stability and hydraulic conductivity based on Three Gorges Dam in China.

The inverse research based on monitoring data and numerical simulation results demonstrates economy and efficiency. In the inversion analysis, optimization algorithms are widely applied for iterative searches over a range of permeability coefficient values. The optimal value is determined by the objective function minimum [28]. Chi [29] constructed an inverse model for the permeability coefficient of a high core rockfill dam based on Particle Swarm Optimization. Wu [30] proposed an improved greedy sampling method based on model reduction technology to rapidly estimate the hydraulic conductivity of a steady state problem. In recent years, with the rapid development of machine learning, many scholars have attempted to use support vector machines [31] and back propagation artificial neural networks [32] for applications in inverse fields, while a BP network is a multi-layer feedforward network trained by the error backpropagation algorithm and is also one of the most widely used network models. Its learning principle is using the steepest descent method to adjust thresholds and weights to make the network’s sum of squares of errors the minimum. It has been proven that BP networks have great accuracy in predicting parameters on the basis of observed value in various engineering research areas [33,34,35]. The marine predator algorithm is a nature-inspired metaheuristic, applying Lévy and Brownian movements in design, and is confirmed to be a high-performance optimizer [36]. Due to its excellent performance, it can be used to search for the optimal solution during the inversion process to minimize the objective function.

The main contributions of this study are as follows: (1) A method for seepage parameters inversion based on an advanced optimization named the Marine predator algorithm is developed; (2) Three different methods for seepage parameters inversion are compared, and the advantage of the Marine predator algorithm in seepage parameters inversion is proposed.

2. Saturated–Unsaturated Seepage Control Equation

It has been proven that Darcy’s law is not only suitable for saturated soils but also for unsaturated soils. Compared to saturated soils, the seepage parameters for unsaturated soils are not constant but a function of soil saturation. Darcy’s law of unsaturated soils is as follows:

$$v=k(\theta )J$$

(1)

$$v=\frac{K(\theta )\gamma }{\mu }J=K\frac{K_r(\theta )\gamma }{\mu }J$$

(2)

where $k(\theta )$ is the permeability coefficient, $K(\theta )$ is the penetration of unsaturated soils, which is a function of the volumetric moisture content (the volume of water per unit volume), $K$ is the penetration of saturated soils, $K_r(\theta )$ is relative penetration, and $\mu$ is the dynamic viscosity of water.

The two-dimensional saturated–unsaturated control equation is as follows:

$$-\frac{\partial (\rho v_i)}{\partial x_i}+Q=\frac{\partial (\rho \theta )}{\partial t}$$

(3)

where $\rho$ is the density of water and $Q$ is the flow.

According to Darcy’s law, Equation (3) can be rewritten as

$$\frac{\partial }{\partial x_i}[k_{ij}^s{k}_r(h_c)\frac{\partial h_c}{\partial x_j}+k_{i3}^s{k}_r(h_c)]+Q=[C(h_c)+\beta S_s]\frac{\partial h_c}{\partial t}$$

(4)

where $k_{ij}^s$ is the saturated permeability tensor, $k_r$ is the relative permeability coefficient in a saturated area, $k_r$ = 1, $C(h_c)=\frac{\partial \theta }{\partial h_c}$ is the water capacity, indicating changes in the moisture content as the water head changes. $\beta$ is 0 in an unsaturated area and 1 in a saturated area. $S_s$ is the unit storage factor, it can be seen as 0 in most instances.

In a stable seepage field, Equation (4) can be rewritten as

$$\frac{\partial }{\partial x_i}[k_{ij}^s{k}_r(h_c)\frac{\partial h_c}{\partial x_j}+k_{i3}^s{k}_r(h_c)]+Q=0$$

(5)

There are two common models for obtaining $k_r$, the van Genuchten model and Brooks–Corey model. The van Genuchten model is as follows:

$$k_r({\theta }_w)={{\theta }_w}^{1/2}{[1-{(1-{{\theta }_w}^{\frac{1}{m}})}^m]}^2$$

(6)

where ${\theta }_w=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$ is effective saturation, ${\theta }_r$ is the saturated moisture content, ${\theta }_s$ is the remaining moisture content, and $\theta$ is the moisture content.

$${\theta }_w={[\frac{1}{1+{(\alpha h_c)}^n}]}^m$$

(7)

$$k_r(h_c)=\frac{{\{1-{(\alpha h_c)}^{n-1}{[1+{(\alpha h_c)}^n]}^{-m}\}}^2}{{[1+{(\alpha h_c)}^n]}^{\frac{m}{2}}}$$

(8)

where $m=1-\frac{1}{n}$, $\alpha$, and $m$ are related to the soil type.

The Brooks–Corey model is as follows:

$${\theta }_w={(\frac{P_d}{P_c})}^{-\lambda }$$

(9)

where $P_d$ is the initial pressure, $P_c$ is the pore water pressure, and $\lambda$ is the shape parameter.

3. Inversion Analysis Principle

3.1. BP Network

3.1.1. Forward Propagation

In this step, raw input data can be converted to output data in hidden layers. The theory is shown in Figure 1.

Each input datum is multiplied by the weight and the activation function to obtain the hidden layer data. The sigmoid function is a common activation function, which is as follows:

$$S(x)=\frac{1}{1+e^{-x}}$$

(10)

Figure 2 shows the drawing of the sigmoid function; each datum is within 0 to 1. Use some matrix to represent the content of Figure 1, and define the input matrix X, hidden layer1 matrix H1, hidden layer2 matrix H2, and output matrix Y as

$$X=\left(\begin{array}{ccccc}in1& in2& in3& \cdots & in8\end{array}\right)$$

(11)

$$H1=\left(\begin{array}{ccc}H{1}_1& H{1}_2& H{1}_3\end{array}\right)$$

(12)

$$H2=\left(\begin{array}{ccc}H{2}_1& H{2}_2& H{2}_3\end{array}\right)$$

(13)

$$Y=\left(\begin{array}{cc}Ou1& Ou2\end{array}\right)$$

(14)

Define the weight matrix and activation function matrix as

$$U=\left(\begin{array}{ccc}{u}_{11}& u_{12}& u_{13}\\ u_{21}& u_{22}& u_3\\ ⋮& \ddots & ⋮\\ u_{81}& \cdots & a_{83}\end{array}\right)$$

(15)

$$H_0=sigmoid(H1)$$

(16)

$$W=\left(\begin{array}{cc}{w}_{11}& w_{12}\\ w_{21}& w_{22}\end{array}\right)$$

(17)

$$H_0=sigmoid(H2)$$

(18)

$$V=\left(\begin{array}{cc}{v}_{11}& v_{12}\\ v_{21}& v_{22}\\ v_{31}& v_{32}\end{array}\right)$$

(19)

The mathematical expression of the Bp network is as follows:

$$X\times U\to H1\to H_0\times W\to H2\to H_0\times V=Y$$

(20)

The process from X to Y is the forward propagation of the Bp network.

3.1.2. Back Propagation

The information transferred by back propagation is an error between the calculated data and the real data, which can be defined as the loss function with a minimum value. The smaller the error, the closer the predicted data are to the actual value. The process of continuously adjusting the weight $U$, $W$, $V$ is the back propagation by the gradient descent method. After the calculation, the original weight matrix is added to the gradient and the learning rate to obtain the new weight for the next calculation. The theory of the gradient descent method is shown in Figure 3.

The network is updated after a forward propagation and a back propagation, and the training of the network is a replication of the forward propagation and the back propagation. As training progresses, the accuracy of the predictions will become better and better. It can be used to improve the accuracy of forecasts to improve the number of layers and the number of nodes.

3.2. Marine Predator Algorithm

The marine predator algorithm (MPA) is a novel and real nature-inspired optimization method based on the marine predator’s complementary time strategy and optimal encounter problem. It has been shown that many animals and marine creatures follow a Lévy flight pattern as their optimal foraging policy [36,37,38,39,40]. Researchers demonstrate that Lévy evolved as an optimal search policy among predators in response to patch prey distribution. Marine predators use the Lévy strategy for environments with a low concentration of prey, while they employ Brownian movement for the areas with abundant prey. The predator’s strategy is to maximize the probability of encountering prey; from a natural perspective, there are three relationships between predators and prey: (1) the predator is faster than the prey; (2) the predator is slower than the prey; (3) the predator and the prey have almost equal speeds. In each phase, the predator has to use the best movement strategy to reach the prey in the best step size, so the design of the MPA matches the rules of the marine predator’s foraging strategy, and the mathematical model of the MPA is very similar to the ecological model of nature [33].

3.2.1. Mathematical Expression

MPA is a population-based method; the initial solution is uniformly distributed over the search space:

$$X_0=X_{\min }+rand(X_{\max }-X_{\min })$$

(21)

where $X_{\min }$ and $X_{\max }$ are the lower and higher boundary of the variables and $rand$ is a uniform random vector from 0 to 1. Build an $Elite$ matrix to simulate the top predator:

$$Elite={\left(\begin{array}{ccc}{X}_{1,1}^{\mathrm{{\rm I}}}& \dots & X_{1,d}^{\mathrm{{\rm I}}}\\ ⋮& \ddots & ⋮\\ X_{n,1}^{\mathrm{{\rm I}}}& \cdots & X_{n,d}^{\mathrm{{\rm I}}}\end{array}\right)}_{n\times d}$$

(22)

$X^{\mathrm{{\rm I}}}$ represents the top predator vector, which is replicated n times to construct the $Elite$ matrix. After each iteration, update the $Elite$ matrix if there exists a better predator. Another matrix with the same dimension as the $Elite$ matrix is named $\mathrm{Pr}ey$, where the predator updates its position based on it. The $\mathrm{Pr}ey$ matrix is as follows:

$$\mathrm{Pr}ey={\left(\begin{array}{ccc}{X}_{1,1}& \dots & X_{1,d}\\ ⋮& \ddots & ⋮\\ X_{n,1}& \cdots & X_{n,d}\end{array}\right)}_{n\times d}$$

(23)

The process is divided into three phases:

Phase 1: When the predator is moving faster than the prey, the best strategy for the predator is not moving, and the mathematical model is as follows:

$$\begin{array}{l}\overrightarrow{stepsiz{e}_i}=\overrightarrow{R_B}\otimes (\overrightarrow{Elit{e}_i}-\overrightarrow{R_B}\otimes \overrightarrow{\mathrm{Pr}e{y}_i})i=1,\dots ,n\\ \overrightarrow{\mathrm{Pr}e{y}_i}=\overrightarrow{\mathrm{Pr}e{y}_i}+P.\overrightarrow{R}\otimes \overrightarrow{stepsiz{e}_i}\end{array}$$

(24)

where $\overrightarrow{R_B}$ is a vector based on a normal distribution representing the Brownian motion. The multiplication of $\overrightarrow{R_B}$ by $\overrightarrow{\mathrm{Pr}e{y}_i}$ simulates the movement of prey. $P=0.5$ is a constant number, and $\overrightarrow{R}$ is a vector of uniform random numbers in [0, 1]. This scenario happens in the first one-third of iterations when the velocity is high with a high exploration ability. $\mathrm{I}ter$ is the current iteration, while $Max\_Iter$ is the maximum.

Phase 2: In this phase, the velocity of the predator is almost equal to that of the prey. The predator and prey are looking for their prey. This section occurs in the middle of optimization; half of the population is designed to explore, while the other half is to develop. The prey are responsible for development and predators are responsible for exploration. Based on this rule, if the prey moves in Lévy, the best strategy of the predator is Brownian.

$$\frac{1}{3}Max\_Iter<\mathrm{I}ter<\frac{2}{3}Max\_Iter$$

For the first half of the population,

$$\begin{array}{l}\overrightarrow{stepsiz{e}_i}=\overrightarrow{R_L}\otimes (\overrightarrow{Elit{e}_i}-\overrightarrow{R_L}\otimes \overrightarrow{\mathrm{Pr}e{y}_i})i=1,\dots ,n\\ \overrightarrow{\mathrm{Pr}e{y}_i}=\overrightarrow{\mathrm{Pr}e{y}_i}+P.\overrightarrow{R}\otimes \overrightarrow{stepsiz{e}_i}\end{array}$$

(25)

where $R_L$ is a vector of random numbers based on Lévy distribution representing Lévy movement. Since most of the Lévy distribution step size is small, it is helpful for development. For the second half of the population, assume that

$$\begin{array}{l}\overrightarrow{stepsiz{e}_i}=\overrightarrow{R_B}\otimes (\overrightarrow{R_B}\otimes \overrightarrow{Elit{e}_i}-\overrightarrow{\mathrm{Pr}e{y}_i})i=1,\dots ,n\\ \overrightarrow{\mathrm{Pr}e{y}_i}=\overrightarrow{Elit{e}_i}+P.CF\otimes \overrightarrow{stepsiz{e}_i}\end{array}$$

(26)

where $CF={(1-\frac{Iter}{Max\_Iter})}^{(2\frac{Iter}{Max\_Iter})}$ is considered an adaptive parameter for controlling the step size for predator movement. The multiplication of $\overrightarrow{R_B}$ and $Elite$ simulates the movement of the predator in a Brownian manner, while the prey updates its position based on the movement of the predators.

Phase 3: In this phase, the velocity of the predator is faster than that of the prey. This scenario happens in the last phase of the optimization with a high exploration capability. The best strategy for the predator is Lévy movement. It is presented as

While $\mathrm{I}ter>Max\_Iter$,

$$\begin{array}{l}\overrightarrow{stepsiz{e}_i}=\overrightarrow{R_L}\otimes (\overrightarrow{R_L}\otimes \overrightarrow{Elit{e}_i}-\overrightarrow{\mathrm{Pr}e{y}_i})i=1,\dots ,n\\ \overrightarrow{\mathrm{Pr}e{y}_i}=\overrightarrow{Elit{e}_i}+P.CF\otimes \overrightarrow{stepsiz{e}_i}\end{array}$$

(27)

Equation (27) simulates the movement of the predator in the Lévy strategy to help update the prey position.

3.2.2. Eddy Formation and FADs’ Effect

There exists an environmental problem which will change the marine predator’s behavior, such as eddy formation or Fish Aggregating Devices (FADs). Sometimes, a shark will take a longer jump in different dimensions, probably to find an environment with another prey distribution. An FAD is seen as the local optima; in the process of searching, consider taking a longer jump to avoid the stop in the local optima. Thus, FADs’ effect is mathematically presented as follows:

$$\begin{array}{l}\overrightarrow{\mathrm{Pr}e{y}_i}\\ =\left\{\begin{array}{ll}\overrightarrow{\mathrm{Pr}e{y}_i}+CF[\overrightarrow{X_{\min }}+\overrightarrow{R}\otimes (\overrightarrow{X_{\max }}-\overrightarrow{X_{\min }})]\otimes \overrightarrow{U}& \mathrm{if} \mathrm{r}<=\mathrm{FADs}\\ \overrightarrow{\mathrm{Pr}e{y}_i}+[FADs(1-r)+r](\overrightarrow{\mathrm{Pr}e{y}_{r1}}-\overrightarrow{\mathrm{Pr}e{y}_{r2}})& \mathrm{if} \mathrm{r}>=\mathrm{FADs}\end{array}\right\}\end{array}$$

(28)

where $FADs=0.2$ is the probability of FADs’ effect on the optimization process and $\overrightarrow{U}$ is the binary vector with arrays including 0 and 1. This is constructed by generating a random vector in [0, 1] and changing its array to 0 if the array is less than 0.2 and 1 if it is greater than 0.2. $\mathrm{r}$ is the uniform random number in [0, 1]. $\overrightarrow{X_{\min }}$ and $\overrightarrow{X_{\max }}$ are the vectors containing the lower and upper bounds of the dimensions. r1 and r2 subscripts denote random indexes of the prey matrix. Based on this, the marine predator has a good memory and remains in the place where they have successfully foraged. The current iteration solution is compared with an equal solution in a prior iteration, and it is updated if the current iteration is more fitted. The pseudocode of the MPA is shown in Algorithm 1.

Algorithm 1. The pseudocode of the MPA.

Confirm search scope Initialize the search agents (prey) populations $i$ = 1, …, n While termination criteria are not met, calculate the fitness, construct the elite matrix, and accomplish memory saving If iter < Max_iter/3, update prey based on Equation (24) Else if Max_iter/3 < iter < 2* Max_iter/3, for the first half of the populations, update prey based on Equation (25); for the other half of the populations, update prey based on Equation (26) Else if iter > 2* Max_iter/3, update prey based on Equation (27) End if accomplish memory saving and Elite update, applying FADs effect, and update based on Equation (28) End

4. Inversion of the Penetration Parameters of a High Core Rockfill Dam

4.1. Parameters Inversion Process

In this paper, an inversion method of the unsaturated seepage parameters of a core rockfill dam based on the MPA and a BP neural network is proposed. The flow chart is shown in Figure 4, and the process is described as follows:

Step 1: Construct the typical 2D cross-section of the actual dam as the actual size. Set up dam partitions and set different materials with different calculation parameters. Input the storage water level process data.

Step 2: Mesh the model, perform a numerical simulation of the model, and obtain the results of the set permeability coefficient combination.

Step 3: Calculate a series of pore pressure with different seepage parameters combinations, and set the rock’s as k and the core’s as K.

Step 4: Train the BP neural network by calculating the data as in step 3, set up upper and lower bounds of k and K, use the MPA method to find a solution, make BP artificial neural network predictions, and find the k and K that make the objective function the minimum.

Step 5: Take the observed value into the trained BP neural artificial network, and obtain the predicted value directly.

Step 6: Add the MPA into the MATLAB–COMSOL linked program, and obtain the k and K that make the model’s results’ objective function the minimum.

Step 7: Ass the k and K obtained from Step 4, Step 5, and Step 6 to the finite element model to calculate the pore pressure of the observer, calculate the absolute error and relative error, and compare the difference between three methods.

Step 8: Bring the parameters into the working conditions of different water storage speeds for the calculation to verify the accuracy of the calculation of the parameters obtained.

4.2. An Actual Dam for Finite Element Analysis

A hydropower station is located at Sichuan Province in China, which is one of the key projects for cascade hydropower development in the Daduhe river basin. This reservoir is the controlling reservoir of the upper reaches of the main stream. The normal water level is 2500.00 m, the design flood level is 2501.81 m, the check flood level is 2504.96 m, and the capacity under a normal water level is 2732 million m³. The dam type is a high core rockfill dam with a maximum dam height of 314 m, a span of 1300 m along the river, and a span of 648.66 m along the dam axis. Since it is such a huge project, it is not feasible to take soil from the dam body to conduct permeability coefficient experiments. There are many pore water pressure detectors buried in the dam body. Based on the pore pressure detector data, the permeability parameters can be inverted to judge the operation of the dam. The material division, the meshing, and the distribution of the observer are as follows.

The model in Figure 5 of the actual dam and the sand layer is built as its actual size, the upstream and downstream ranges are each taken as 0.5 times the dam height, and the bedrock depth is twice the dam height. The calculation of the dam body pore water pressure is carried out in COMSOL 6.0 FEM software, which can accurately calculate the trend of the wetting line inside the dam body and the distribution of pore water pressure in the dam body. The meshing is shown in Figure 6. The distribution of measuring points in the core wall is shown in Figure 7.

4.3. The Initial Calculated Seepage Parameter

It has been shown that Darcy’s law is also adapted in unsaturated flow; in the water storage stage, some parts of the dam body are still in an unsaturated state, so it is necessary to carry out saturated–unsaturated seepage analysis. It can be seen in Equation (8) that in the van Genuchten model (hereinafter referred to as the vg model), the permeability coefficient of porous media in the unsaturated state is jointly determined by the water content, the saturated permeability coefficient and the shape parameters $\alpha$ and $n$. The default calculated permeability coefficient is as

Table 1. The default calculated permeability coefficient.

	k(cm/s)	$\mathit{\alpha }$	$\mathit{n}$
gravel core wall	7.00 × 10−6	0.02	2
upstream rock	1.00 × 10−1	0.02	2
downstream rock	1.00 × 10−1	0.02	2
curtain grouting	1.00 × 10−7	0.02	2
transition layer	3.00 × 10−2	0.02	2
filter layer	8.00 × 10−3	0.02	2
sand layer	3.00 × 10−3	0.02	2
bed rock	1.00 × 10−5	0.02	2

, $\alpha$ = 0.02 and $n$ = 2:

There are three stages for the initial water storage of the reservoir. The initial water storage speed is 2 m/d and lasts for 30 days, the mid-term storage speed is 2 m/d and lasts for 40 days, the late water storage speed is 1 m/d and lasts for 80 days, the total water storage period is 150 days, and the water storage elevation is 220 m. The water storage process line is as Figure 8:

In the transient analysis, the total time length is 150 days of total water storage time, and the step size is 1 day. In the first calculation with preset penetration parameter values, the observed result and calculated result are as follows (in this paper, several of the observation points are analyzed: 1, 2, 3, 4, 5, 10, 14).

It can be seen in

Table 2. The result of the first calculation with preset penetration parameter values.

Observation Point Number	Pore Water Pressure Observed (m)	Pore Water Pressure Calculated (m)
1	5.81	6.41
2	19.82	18.81
3	30.24	28.84
4	46.63	45.76
5	53.62	59.81
10	161.51	163.85
14	245.35	244.09

that the pore water pressure shows an increasing trend from top to bottom in the dam body. However, there exists a large error between the calculated value and the observed value at each observation point, so it is necessary to conduct inversion to understand the seepage parameters in the dam body. The distribution of the pressure water head is as Figure 9.

It can be seen based on the foregoing that the pressure water head distributed in the dam body is connected to the material’s shape parameters, moisture content, and saturation seepage parameters. Since the case simulates the storage process of a core rockfill dam without the impact of rainfall infiltration, the main influence factors are shape parameters and saturation seepage parameters. Since the upstream rock and gravel core wall are different, it is not accurate enough to allow them to be calculated using the same shape parameters. Figure 1 shows that the filter layer and transition layer have little influence on pressure head changes, and the effects of the two are ignored for now. In enumerating a series of parameters to discuss the impact of shape parameters, the result is as

Table 3. Effects caused by changes in shape materials.

Pressure Water Head (m)	Observe 1	Observe 2	Observe 3	Observe 4	Observe 5	Observe 10	Observe 14
rock a = 0.02 n = 2	6.41	18.81	28.84	45.76	59.81	163.85	244.09
core a = 0.02 n = 2
rock a = 0.02 n = 2	6.6	21	31.55	47.23	58.78	162.07	243.87
core a = 0.03 n = 2
rock a = 0.03 n = 2	9.9	21.21	30.11	47.28	61.76	164.11	244.15
core a = 0.02 n = 2
rock a = 0.03 n = 2	9.87	21.1	28.52	45.01	58.1	162.35	243.95
core a = 0.03 n = 2

:

Compared with Table 2, it can be seen that on the upper part of the core wall, increasing the rock’s shape material $\alpha$ will have a great impact; however, in another area, the impact caused by the rock and core is not much different. It is worth nothing that changing the parameters of both at the same time has little impact on the pressure water head apart from the upper area; in this paper, $\alpha$ = 0.02 and n = 2 are taken for the subsequent inversion analysis of the rock and core wall. Assume the rock’s permeability ranges from 0.06 cm/s to 0.14 cm/s and the core wall’s permeability ranges from 3.00 × 10⁻⁶ cm/s to 11.0 × 10⁻⁶ cm/s, and assume a series of computational combinations for BP neural network training and MPA seek optimization. The maximum and minimum values are as

Table 4. The maximum and minimum values of the pressure water head of computational combinations.

K (cm/s)	K (cm/s)	Observe 1	Observe 2	Observe 3	Observe 4	Observe 5	Observe 10	Observe 14
3.00 × 10−6	0.1	3.75	9.68	12.13	28.53	38.91	156.86	243.22
1.10 × 10−6	0.14	6.56	22.3	33.97	51.98	67.73	165.41	244.31

:

It has been proven that the calculation range includes the true permeability coefficient value of the rock and core wall; this method can accurately predict the permeability coefficient value. Set the objective function as follows:

$$f_n={\displaystyle \sum _i^7{(\frac{c_i-o_i}{{(o_i)}^2})}^{\frac{1}{2}}}$$

(29)

where $f_n$ is the objective function value calculated in each iteration, $n$ is the iteration, the maximum number of iterations is taken as 500, $n_{\max }$ = 500, $c_i$ is the calculated value of one observation point, and $o_i$ is the observed value of one observation point. The combination of permeability coefficients that minimizes the objective function $f_n$ value is the required permeability coefficient. The permeability coefficient calculated by BP–MPA is as follows: K = 4.016 × 10⁻⁶ cm/s, k = 0.06980 cm/s, and the minimum value of the objective function is 0.3438. Substitute the observed water head into the trained BP artificial neural network model, predicting that K = 6.93 × 10⁻⁶ cm/s, k = 0.442 cm/s, and the minimum value of the objective function is 0.3124. Apply the MPA algorithm to COMSOL with the MATLAB program; after 500 iterations, it can be seen that K = 6.89 × 10⁻⁶ cm/s, k = 0.949 cm/s, and the minimum value of the objective function is 0.3493 (K is the permeability of the core; k is the permeability of the rock).

4.4. Verification of the Permeability Coefficient

Add the permeability coefficient into the finite element method to calculate the pressure water head of the observation point to verify the accuracy. The absolute error and relative error are as Figure 10:

It can be seen that the calculation errors of the observation points at the lower part of the core wall are relative to each other. In other areas, the three inversion methods have their own advantages at different computing nodes. The average of the seven-point absolute errors and relative errors is as

Table 5. The error’s average and variance.

	Average of Absolute Error (m)	Average of Relative Error	Variance in Absolute Error	Variance in Relative Error
MPA-BP	3.52	0.105	3.15	0.109
BP	2.21	0.058	1.43	0.035
COMSOL–MATLAB	3.68	0.083	3.78	0.068

:

It can be seen that the BP method’s average of the absolute error is the minimum, though the absolute error of the MPA-BP method is smaller than that of the COMSOL–MATLAB method, the second’s relative error is smaller, and the law of variance is also similar to the law of error. So, the COMSOL–MATLAB method may cause a large deviation at a certain code, and the BP method shows the best accuracy in terms of prediction.

There exists another storage plan that involves adding the calculated permeability coefficient into the second calculated conditions and changing the medium-term storage speed to 1 m/d; the result is as

Table 6. Result of the second calculation condition.

Pressure Water Head (m)	K (cm/s)	K (cm/s)	1	2	3	4	5	10	14
observed value			5.9	17.9	30.4	46.4	58.1	152	247.9
MPA-BP	4.02 × 10−6	0.0698	5.64	13.73	19.84	36.64	49.87	160.98	243.66
BP	6.93 × 10−6	0.442	6.23	19.52	28.98	46.68	61.46	164.56	244.46
COMSOL–MATLAB	6.89 × 10−6	0.949	5.64	19.05	29.47	45.84	60.64	164.56	244.49

:

The objective function calculated by the permeability coefficient obtained by COMSOL–MATLAB is the minimum, which is 0.253. The distribution of the pore water pressure in the dam body is as Figure 11:

Compared with Figure 9, the wetting line in the core wall is closer to the downstream direction, since the extension of the water storage time causes the materials in the dam body to become more saturated.

5. Discussions

The seepage parameters of the dam material will change as the water storage progresses, affecting the safe operation of the dam. Many inversion methods have been proposed for accurately obtaining the dam’s seepage parameters based on the observed value. The BP artificial network has been widely used in hydrological prediction, but it can also be used in the prediction of seepage parameters. The MPA is a nature-inspired metaheuristic, applying Lévy and Brownian movements in its design, and it is confirmed to be a high-performance optimizer, so it can also be used in seepage parameters inversion.

This paper presents an inversion analysis of a dam’s permeability coefficient. It benefits from the great relationship between the calculation software MATLAB 2018 and the FEM software COMSOL 6.0, it takes the MPA into the inversion-forward calculation loop iteration, and it obtains the optimal solution. Compared with the results obtained by the BP artificial neural network and BP–MPA method, the BP method result’s error is the smallest. But in some work conditions, the MPA shows better accuracy in permeability coefficient prediction, since a BP artificial neural network needs massive data for network training, while an MPA can find global optimization based on a small amount of data. So, in many conditions, the MPA is applicable for the inversion of a dam’s seepage parameters. It is hoped that more research regarding the comparison of different algorithms can be conducted in order to perform parameter inversion more efficiently to ensure the safety of dam operation.

6. Conclusions

This paper proposed a new method for the inverse permeability coefficient of a core-rock dam named the marine predator algorithm, which is a metaheuristic optimization algorithm added to a 2D model calculation of a core-rock dam in China. Incorporating the BP artificial neural network into the comparison, network training is conducted by setting up 81 datasets and predicting the permeability coefficient. Three inversion methods are compared: the MPA is added to the BP artificial network to calculate the optimal solution, the BP network is used to predict the real permeability coefficient, and a piece of the MPA MATLAB code is written to link COMSOL 6.0 software performing inverse-forward analysis. The conclusions are as follows:

(1) In the first calculation condition, the BP artificial network algorithm shows better accuracy in terms of prediction results due to the large dataset for training; in some observed points, the other two methods have smaller errors. The mean error and variance of the absolute error of the BP artificial neural network algorithm are both the smallest, which are 2.21 m and 1.43 m, respectively. The mean error and variance of the relative error of the BP artificial neural network algorithm are both the smallest, which are 0.058 and 0.035, respectively. The mean error and variance of the absolute error of BP–MPA and COMSOL–MATLAB are 3.52 m, 3.15 m and 3.68 m, and 3.78 m. The mean error and variance of the relative error of BP–MPA and COMSOL–MATLAB are 0.105, 0.109 and 0.083, and 0.068. This shows that the second method may lead to fewer discrete errors, while the specific error may be abnormal.

(2) There exists another calculation condition that has a slower water storage speed; in this condition, add the permeability coefficient obtained by the MATLAB link COMSOL method into the finite element model. Its error is much smaller than the results calculate by the permeability coefficient and predicted by the BP artificial neural network since the network lacks data for training.

(3) The MPA is reliable and efficient in inversion when there is only a small amount of data. When there are more data, conducting network training may lead to a higher accuracy in predicting inversion parameters.