**Author details**

Hemali Rathnayake\* and Sheeba Dawood Department of Nanoscience, Joint School of Nanoscience and Nanoengineering, University of North Carolina at Greensboro, Greensboro, NC, USA

\*Address all correspondence to<EMAIL_ADDRESS>
© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**161**

*Coordination Polymer Frameworks for Next Generation Optoelectronic Devices*

[10] Liu Y, Zhao Y, Chen X (2019) Bioengineering of metal-organic frameworks for nanomedicine.

[11] Rowsell JL, Yaghi OM (2004) Metal–organic frameworks: a new class of porous materials. Microporous and

[12] Kotzabasaki M, Froudakis GE (2018) Review of computer simulations

[13] Rowsell JL, Spencer EC, Eckert J, Howard JA, Yaghi OM (2005) Gas adsorption sites in a large-pore metal-organic framework. Science

[14] Yang J, Trickett CA, Alahmadi SB, Alshammari AS, Yaghi OM (2017) Calcium L-lactate frameworks as naturally degradable carriers for pesticides. Journal of the American Chemical Society 139:8118-8121

[15] Silva CG, Corma A, García H (2010) Metal–organic frameworks as semiconductors. Journal of Materials

[16] Stassen I, Burtch N, Talin A,

Falcaro P, Allendorf M, Ameloot R (2017) An updated roadmap for the integration of metal–organic frameworks with electronic devices and chemical sensors. Chemical Society Reviews 46:3185-3241

[17] Usman M, Mendiratta S, Lu K-L (2017) Semiconductor metal–organic frameworks: Future low-bandgap materials. Advanced Materials

[18] Stavila V, Talin AA, Allendorf MD (2014) MOF-based electronic and optoelectronic devices. Chemical Society

Chemistry 20:3141-3156

29:1605071

Reviews 43:5994-6010

on anti-cancer drug delivery in MOFs. Inorganic Chemistry Frontiers

mesoporous materials 73:3-14

Theranostics 9:3122

5:1255-1272

309:1350-1354

*DOI: http://dx.doi.org/10.5772/intechopen.94335*

[1] Zafar F, Sharmin E (2016) Metal-Organic Frameworks. BoD–Books on

[2] Ogawa T, Iyoki K, Fukushima T, Kajikawa Y (2017) Landscape of research areas for zeolites and metal-organic frameworks using computational classification based on citation networks. Materials 10:1428

[3] Furukawa H, Cordova KE, O'Keeffe M, Yaghi OM (2013) The chemistry and applications of metalorganic frameworks. Science 341:

Materials 25:27-33

slovenica 53:117

[4] Eslava S, Zhang L, Esconjauregui S, Yang J, Vanstreels K, Baklanov MR, Saiz E (2013) Metal-organic framework ZIF-8 films as low-κ dielectrics in microelectronics. Chemistry of

[5] Yap MH, Fow KL, Chen GZ (2017) Synthesis and applications of MOFderived porous nanostructures. Green Energy & Environment 2:218-245

[7] Espallargas GM, Coronado E (2018) Magnetic functionalities in MOFs: from the framework to the pore. Chemical

[6] Logar NZ, Kaucic V (2006) Nanoporous materials: from catalysis and hydrogen storage to wastewater treatment. Acta chimica

Society Reviews 47:533-557

[8] Vikrant K, Kumar V, Kim K-H, Kukkar D (2017) Metal–organic frameworks (MOFs): potential and challenges for capture and abatement of ammonia. Journal of Materials Chemistry A 5:22877-22896

[9] Ma M, Lu L, Li H, Xiong Y, Dong F (2019) Functional metal organic framework/sio2 nanocomposites: From versatile synthesis to advanced applications. Polymers 11:1823

**References**

Demand

*Coordination Polymer Frameworks for Next Generation Optoelectronic Devices DOI: http://dx.doi.org/10.5772/intechopen.94335*


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

**Yongkang Shu 1, Zhenzhong Shen 1,2, Liqun Xu 1,\*, Junrong Duan 1, Luyi Ju 1 and Qi Liu 3**

**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.

**Keywords:** inverse analysis; hydraulic conductivities; Gray Wolf Optimizer

#### **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–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

**Citation:** Shu, Y.; Shen, Z.; Xu, L.; Duan, J.; Ju, L.; Liu, Q. Inverse Modeling of Seepage Parameters Based on an Improved Gray Wolf Optimizer. *Appl. Sci.* **2022**, *12*, 8519. https://doi.org/10.3390/ app12178519

Academic Editors: Danial Jahed Armaghani, Yixia Zhang, Pijush Samui, Ahmed Hussein Kamel Ahmed Elshafie and Aydin Azizi

Received: 14 July 2022 Accepted: 24 August 2022 Published: 25 August 2022

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–12] and ant colony optimization [13–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–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 metaheuristics. 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 *α*, *β*, and *γ* represent the three dominant wolves in rank order. The mathematical model of gray wolf hunting is established as described in Equation (1).

$$\begin{cases} \,d = \left| \mathbf{C} \cdot \mathbf{X}\_p(t) - \mathbf{X}(t) \right|\\ \mathbf{X}(t+1) = \mathbf{X}(t) - \mathbf{A} \cdot \mathbf{d} \end{cases} \tag{1}$$

where *d* is the perceived distance between the gray wolf and the prey; *t* represents the number of iterations. *C* denotes the coefficient vector; *Xp*(*t*) and *X*(*t*) are the positions of the wolf and the prey, respectively; *X*(*t* + 1) 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).

$$\begin{cases} A = a \cdot (2m\_1 - 1) \\ \mathbb{C} = 2m\_2 \end{cases} \tag{2}$$

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).

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

where *X*(*t* + 1) denotes the position of the wolf *ϕ* after update; *Xi*(*i* = 1, 2, 3) represent the position vectors of the wolf *α*, *β*, and *γ*, respectively.

The position vectors of *α*, *β*, and *γ* are expressed by Equation (4).

$$\begin{cases} X\_1 = |X\_a - A\_1 \cdot d\_1| \\ X\_2 = |X\_\beta - A\_2 \cdot d\_2| \\ X\_3 = |X\_\gamma - A\_3 \cdot d\_3| \end{cases} \tag{4}$$

where *Xi*(*i* = 1, 2, 3) stand for the position vectors of the wolf *α*, *β*, and *γ*, respectively; *Xj*(*j* = *α*, *β*, *γ*)represent the prey position perceived by *α*, *β*, and *γ*, respectively. *Ai*(*i* = 1, 2, 3) are the coefficient vectors; *di*(*i* = 1, 2, 3) mean the distance between the three dominant wolves and the prey.

*di*(*i* = 1, 2, 3) could be expressed as Equation (5) follows.

$$\begin{cases} d\_1 = \left| \mathbf{C}\_1 \cdot \mathbf{X}\_{\mathcal{A}} - \mathbf{X}(t) \right| \\ d\_2 = \left| \mathbf{C}\_2 \cdot \mathbf{X}\_{\mathcal{S}} - \mathbf{X}(t) \right| \\ d\_3 = \left| \mathbf{C}\_3 \cdot \mathbf{X}\_{\mathcal{I}} - \mathbf{X}(t) \right| \end{cases} \tag{5}$$

where *Ci*(*i* = 1, 2, 3) are the coefficient vectors; *X*(*t*) 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).

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

where Δ*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).

$$R\_j = lb + \frac{ub - lb}{S}j + rand() \cdot \frac{ub - lb}{S}(j = 1, 2, \dots, S) \tag{7}$$

where *Rj* 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 |*A*| < 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 |*A*| > 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).

$$a(t) = r\_1(\frac{t}{T})^3 + r\_2(\frac{t}{T})^2 + r\_3\frac{t}{T} + r\_4\tag{8}$$

where *a* is the convergence factor; *ri*(*i* = 1, 2, 3, 4) 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).

$$\begin{cases} a(0) = 2, a(T) = 0, a\left(\frac{T}{2}\right) = 1\\ \quad a'(0) < 0, a'(T) < 0 \end{cases} \tag{9}$$

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

$$a(t) = \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$$

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 × 10−3, −1 × 10−3, and −2 × 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 × 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 × 10−3 was chosen to ensure the property of jumping out of the local optimum in the study.

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

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).

$$\begin{cases} \rho\_1 = \frac{\|X\_1\|}{\|X\_1\| + \|X\_2\| + \|X\_3\|}\\ \rho\_2 = \frac{\|X\_2\|}{\|X\_1\| + \|X\_2\| + \|X\_3\|}\\ \rho\_3 = \frac{\|X\_1\| + \|X\_2\| + \|X\_3\|}{\|X\_1\| + \|X\_2\| + \|X\_3\|} \end{cases} \tag{11}$$

Here *ρi* (*i* = 1, 2, 3) are the weighting factors of the first-three level wolves, respectively; *Xi* (*i* = 1, 2, 3) are the Euclidean norms of the first three levels, respectively.

The weights of the wolf *α*, *β*, and *γ* 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.

$$X(t+1) = \frac{6\rho\_1 \cdot X\_1 + 3\rho\_2 \cdot X\_2 + \rho\_3 \cdot X\_3}{10} \tag{12}$$

where *ρi* (*i* = 1, 2, 3) are the weighting factors of the first-three level wolves, respectively; *Xi* (*i* = 1, 2 , 3) represent the position vectors of the wolf *α*, *β*, and *γ* 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.

**Table 2.** Results of the numerical experiments.

**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.

#### **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 *K* = [*k*1, *k*2,..., *kn*] the combination of hydraulic conductivity, in which *ki* represents the hydraulic conductivity of the *i*th media. The number of piezometers and measuring weirs are indicated by *m* and *n*, respectively. *H* = [*H*1, *H*2,..., *HM*] is expressed as a sequence of hydraulic head measurements. Similarly, *Q* = [*Q*1, *Q*2,..., *QN*] is a series of leakage volume measurements. *Hi*(*K*) and *Qj*(*K*) 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).

$$\min f = \left(\sum\_{i=1}^{M} \frac{\|H\_i(\mathbf{K}) - H\_i\|\_2^2}{\|H\_i\|\_2^2}\right)^{\frac{1}{2}} + w \left(\sum\_{i=1}^{N} \frac{\|Q\_j(\mathbf{K}) - Q\_j\|\_2^2}{\|Q\_j\|\_2^2}\right)^{\frac{1}{2}} \tag{13}$$

$$\text{s.t. } \mathbf{K}\_{\text{min}} \le \mathbf{K} \le \mathbf{K}\_{\text{max}}$$

Here *Kmin* and *Kmax* 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, *Kmin* and *Kmax*, 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 *α*, *β*, and *γ*.

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.

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

#### **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.

**Figure 4.** The location of the Kakiwa Dam.

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.

**Figure 5.** The max 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 *PDB*−13 is located downstream of the curtain. *PDB*−24 − *PDB*−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.

**Table 3.** Hydraulic conductivity of stationary medium.

#### *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 *PDB*−24 is approximately 2776 m, and the difference among the values of piezometers *PDB*−24 to *PDB*−27 is not significant, indicating the efficiency of the impermeable curtain.

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

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 m3/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.

**Figure 7.** 2D FE mesh of Kakiwa dam.

#### *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.

**Table 4.** Hydraulic conductivity of each medium determined by IGWO.

#### 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):

$$
\delta\_H = \frac{|H\_i - H|}{\Delta H} \times 100\% \tag{14}
$$

where *δH* is the relative error of the hydraulic head; *Hi* and *H* represent the simulated and measured values of the hydraulic head, respectively. Δ*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 *PDB*−24 − *PDB*−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.

**Table 5.** Comparison between the calculated and measured hydraulic head values.

#### 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).

$$Q\_i = \iint vB dx dz \tag{15}$$

where *Qi* 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).

$$
\delta\_Q = \frac{|Q\_i - Q|}{Q} \times 100\% \tag{16}
$$

where *δQ* is the relative error of the dam body leakage; *Qi* 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.

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

#### **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 m3/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.

**Funding:** This research was funded by the "National Key R&D Program of China, grant number 2019YFC1510802".

**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.

#### **References**


*2.4. Kinematic Wave*

For most simple river systems, the kinematic wave formulation is applicable. The kinematic wave calculates flow wave dissemination and flow velocity variations over a network of streams. The roughness of the slope and river bottom affect kinematic flow waves. By changing the Manning equation (Equation (1)) in the continuity equation and differentiating cross-sectional area with respect to time, the kinematic wave differential equation may be produced (Equation (2)).

$$Q = \frac{1}{n} \frac{A^{5/3}}{B^{2/3}} S\_0^{1/2} \tag{1}$$

$$\frac{\partial Q}{\partial \mathbf{x}} + \alpha \beta Q^{\beta - 1} \frac{\partial Q}{\partial t} = 0 \tag{2}$$

where *Q* is the flow rate, *x* is longitudinal distance in channel, *t* is time, *α* and *β* are functions of hydraulic coefficients, *S*0 is the slope of the channel, *n* is the Manning friction factor, *A* is the cross-sectional area, and *B* is the channel width. The depth exponent is used along with segment geometries to estimate hydraulic coefficients, which are later used to calculate segment flow depths under specific flow rates. [38] gave a complete description of WASP's stream transport. To quantify variations in velocities, widths, and depths across the network, the kinematic wave based on solutions of one-dimensional continuity equations, and a condensed version of the momentum equation that takes the effects of friction and gravity into account, was utilized.

#### *2.5. Discretization*

The river network was discretized using a 500 m interval. This led to 319 horizontal, one-dimensional water segments in WASP. Segment widths, average depths, and depth exponents were among the hydrodynamic metrics produced by HEC-RAS. The depth exponent controls channel shape in WASP. For these segments, a depth exponent value of 0.3 was used to simulate an uneven cross-section. Along with segment widths, slopes, and roughness factors, these depths are used in simulations to determine segment depths. The model's input variables include channel geometry, flow routing, boundary conditions, environmental time functions, loads, and initial segment conditions. Hydraulic parameters were obtained using a calibrated and validated HEC-RAS model and utilized as inputs for the flow functions in WASP [38]. From the 319 HEC-RAS cross-sections separated by 500 m, geometry for the WASP segments (average depths, widths, and slopes) was derived. Based on the station's water quality data closest to the study region, spatial linear interpolation was used to identify the boundary and initial conditions for the current investigation.

The study's duration was from 2020 to 2021. Initial heavy metal concentrations for each model component were derived by interpolating between two nearby long-term monitoring stations: the upstream boundary conditions were defined using the nearest long-term heavy metal data, which was Outlook. The boundary data for downstream was used based on the data from Weldon Ferry.

A segment length of 500 m was used for modelling. This length is deemed appropriate based on processing times, a reasonable uniform volume for segments, and an acceptable mixing. Benthic segments were inserted underneath each section of surface water to quantify erosion and sedimentation. The measured concentration at the upstream and downstream boundaries was used to calibrate the deposition rate of metals. In the HEC-RAS model, the spaces between each pair of cross-sections were split into individual WASP segments. Output files of the hydrodynamic model (HEC-RAS) were specific for each crosssection. These results then needed to be averaged between two consecutive cross-sections to obtain the exact values for each WASP segment. This procedure is shown in Figure 2 and was implemented in Microsoft Excel (Version 2211) for all input values. *Water* **2023**, *15*, x FOR PEER REVIEW 7 of 16

**Figure 2.** Model segmentation. **Figure 2.** Model segmentation.

session.

#### *2.6. Visualization, Calibration, and Validation 2.6. Visualization, Calibration, and Validation*

The modelled heavy metal concentrations and hydrodynamics were calibrated and validated. The calibration of heavy metal concentrations involved several iterations in the parameter space within reasonable parameters of the observed values. The best curve fitting with the observed data was achieved by choosing the optimum parameter values. To validate the model, a run with these values was compared against an independent set of field data. The modelled heavy metal concentrations and hydrodynamics were calibrated and validated. The calibration of heavy metal concentrations involved several iterations in the parameter space within reasonable parameters of the observed values. The best curve fitting with the observed data was achieved by choosing the optimum parameter values. To validate the model, a run with these values was compared against an independent set of field data.

how many geographical grids, x/y plots, or even model result files a user can employ in a

An effective way to review model simulations and calibrate them using collected data is to utilize the post-processor (MOVEM). Results from every WASP run, as well as others, can be visualized using MOVEM. MOVEM allows the modeller the choice of two graph-

An effective way to review model simulations and calibrate them using collected data is to utilize the post-processor (MOVEM). Results from every WASP run, as well as others, can be visualized using MOVEM. MOVEM allows the modeller the choice of two graphical representations of the results: geometric grid and x/y plots. There is no restriction on how many geographical grids, x/y plots, or even model result files a user can employ in a session. *Water* **2023**, *15*, x FOR PEER REVIEW 8 of 16

#### **3. Results 3. Results**

#### *3.1. Deposition Rates and Overall Concentration Trends 3.1. Deposition Rates and Overall Concentration Trends*

First, it was assumed that the deposition rate for copper and nickel should be the same as sediment. Deposition rates for copper and nickel were calibrated using the sampled concentration at the upstream and downstream model boundaries (Figure 3). A previously developed model for sediment transport was used, and the model was calibrated with an acceptable range of deposition [38]. Concentration plots (Figure 4) show the trend for both trace metals at selected sampling sites (Downtown, Clarkboro Ferry, and Highway 312) for both water and sediment samples. Concentration values for all sites can be found in the Supplementary Information. First, it was assumed that the deposition rate for copper and nickel should be the same as sediment. Deposition rates for copper and nickel were calibrated using the sampled concentration at the upstream and downstream model boundaries (Figure 3). A previously developed model for sediment transport was used, and the model was calibrated with an acceptable range of deposition **[**38**]**. Concentration plots (Figure 4) show the trend for both trace metals at selected sampling sites (Downtown, Clarkboro Ferry, and Highway 312) for both water and sediment samples. Concentration values for all sites can be found in the Supplementary Information.

**Figure 3.** Calibrated deposition curve for both trace metals Cu and Ni. **Figure 3.** Calibrated deposition curve for both trace metals Cu and Ni.

Both trace metals were modelled for the same time periods for water and sediment in the years 2020 and 2021 and matched to the field sampling of water and sediment samples to maintain uniformity.

#### *3.2. Copper*

In 2020, copper concentrations were observed in the range of 1–5 µg/L except for two sites (Highway 3 and Weldon Ferry), with values at 42.2 µg/L and 20.2 µg/L, respectively (refer to Supplementary Information). Cu concentrations in river water in 2021 were in the range of 0.8 µg/L to 5 µg/L. That the observed concentrations remained uniform, for the most part, is indicative of the complex mixing of copper across the sampling sites, whereas elevated concentrations of copper at Highway 3 and Weldon Ferry could be attributed to various constructions around the location (bridges and roadway) in this area (S. Prajapati, pers. obs.). Many of the measured concentrations were over the long-term CCME water quality guideline for Cu of 4 µg/L at the given water hardness > 180 mg/L CaCO3, which indicates that risks to aquatic life could not be excluded. Sediment Cu concentrations ranged from <LOQ to 19 mg/kg in 2020 and <LOQ to 15 mg/kg in 2021 [27]. Since the

resulting concentrations are below the CCME ISQG of 35.7 mg/kg, the authors do not expect adverse effects. In sediments, Cu ranged from 1–15 mg/kg, with Clarkboro Ferry being on the higher end for both sampling years, with values in the range of 12–15 mg/kg. *Water* **2023**, *15*, x FOR PEER REVIEW 9 of 16

**Figure 4.** Concentrations of Cu and Ni in water and sediment for selected sites during the 2020 and 2021 field seasons. **Figure 4.** Concentrations of Cu and Ni in water and sediment for selected sites during the 2020 and 2021 field seasons.

Both trace metals were modelled for the same time periods for water and sediment in the years 2020 and 2021 and matched to the field sampling of water and sediment samples to maintain uniformity. *3.2. Copper* In 2020, copper concentrations were observed in the range of 1–5 µg/L except for two Figures 5 and 6 show a longitudinal comparison of measured and modelled copper concentrations in water and sediment for 2020 and 2021. The concentrations modelled through the sampling time points almost uniformly increased between the upstream and downstream boundaries of the study area. Copper concentrations for water samples remained mostly consistent for June 2020, August 2020, and October 2020, but increased gradually and became more consistent in June 2021 and August 2021. Within the study

sites (Highway 3 and Weldon Ferry), with values at 42.2 µg/L and 20.2 µg/L, respectively (refer to Supplementary Information). Cu concentrations in river water in 2021 were in

period, discharge and flow velocity of the South Saskatchewan River were highest in June 2020, which can explain this observed pattern of copper concentrations, as increased discharge would cause significant erosion of riverbanks and stream beds. Thereafter, longitudinal profiles for August 2020, October 2020, and June 2021 remained relatively consistent throughout the segments. *Water* **2023**, *15*, x FOR PEER REVIEW 11 of 16

**Figure 5.** Longitudinal profiles of water concentrations with observed values for both trace metals in the 2020 and 2021 field seasons. **Figure 5.** Longitudinal profiles of water concentrations with observed values for both trace metals in the 2020 and 2021 field seasons.

**Figure 6.** Longitudinal profiles of sediment concentrations for both trace metals in the 2020 and 2021 field seasons. **Figure 6.** Longitudinal profiles of sediment concentrations for both trace metals in the 2020 and 2021 field seasons.

#### *3.3. Nickel*

*3.3. Nickel* In 2020, nickel concentrations ranged consistently from 2–6 µg/L for water samples and ranged from 0.3–1.4 µg/L in 2021. The high-flow season in 2020 (due to elevated discharge and precipitation mentioned above) helps explain the almost three- to fourfold greater nickel concentrations in that year. Levels of nickel in sediments for 2020 ranged from 10.5–25 mg/kg and mostly remained consistent throughout the year. However, in contrast to concentration values for water samples, the concentration of nickel in 2021 remained around the same range, with a very slight decrease, ranging from 10–20 mg/kg. No measured values were observed to be over the long-term CCME water quality guide-In 2020, nickel concentrations ranged consistently from 2–6 µg/L for water samples and ranged from 0.3–1.4 µg/L in 2021. The high-flow season in 2020 (due to elevated discharge and precipitation mentioned above) helps explain the almost three- to fourfold greater nickel concentrations in that year. Levels of nickel in sediments for 2020 ranged from 10.5–25 mg/kg and mostly remained consistent throughout the year. However, in contrast to concentration values for water samples, the concentration of nickel in 2021 remained around the same range, with a very slight decrease, ranging from 10–20 mg/kg. No measured values were observed to be over the long-term CCME water quality guidelines for Ni of 150 µg/L at the given water hardness > 180 mg/L CaCO3.

lines for Ni of 150 µg/L at the given water hardness > 180 mg/L CaCO3. Longitudinal profiles for nickel (Figures 5 and 6) followed almost the same trend as those for copper, with very little variability in concentration values during the different Longitudinal profiles for nickel (Figures 5 and 6) followed almost the same trend as those for copper, with very little variability in concentration values during the different sampling time points. Nickel concentrations remained consistent for August 2020, October 2020, June 2021, and August 2021, with a slight gradual increment with every major segment (i.e., sampling station). The concentration values followed a similar trend of gradual increase. For sediment samples, the modelled progressions again followed similar patterns in August 2020, October 2020, and June 2021, with a consistent and slight decline throughout the study area. Again, high flow in 2020 could explain the three- to fourfold greater nickel concentrations in 2020 compared to 2021 water samples. Similarly, for sediment samples, settling and transverse mixing can contribute to consistent and slight declines in the concentration values.


**1. Introduction**

Since the end of the 20th century, environmental concerns and the development of sustainable strategies have been on the rise [1]. This concern carries over to the present day, where plastic represents a large part of human life, with up to 368 million tonnes of plastic produced globally in 2019 [2]. The problem arises mainly because of the negative impact on the ecosystems generated by the amount of plastic waste, since only 173 million tonnes are collected for recycling or landfill [3], but also by the exploitation of non-renewable natural resources, such as petroleum. The problem is magnified in the packaging industry, as it is the largest consumer of plastics, accounting for 39.6% of total plastic consumption in Europe [2]. In view of this situation, research is mainly focused on obtaining polymeric materials with formulations based on renewable resources or even with the property of being biodegradable. A wide variety of biopolymers or bio-based materials have now been obtained that contribute to this sustainable development.

Bio-based polymers are polymers that are made from biological substances, i.e., nonfossil materials [4]. These materials may or may not be biodegradable depending on their ability to be broken down by microorganisms [4]. Currently, some of the most promising biopolymers are extracted from biomass (starch, cellulose, protein, chitin, etc.), such as thermoplastic starches (TPSs), given the abundance of polysaccharides in the biomass. Another promising alternative are also those obtained by microbial production, such

#### **Citation:** Perez-Nakai, A.;

Lerma-Canto, A.; Domingez-Candela, I.; Garcia-Garcia, D.; Ferri, J.M.; Fombuena, V. Comparative Study of the Properties of Plasticized Polylactic Acid with Maleinized Hemp Seed Oil and a Novel Maleinized Brazil Nut Seed Oil. *Polymers* **2021**, *13*, 2376. https://doi.org/10.3390/ polym13142376

Academic Editors: Andreia F. Sousa, José Miguel Ferri, Vicent Fombuena Borràs and Miguel Fernando Aldás Carrasco

Received: 30 June 2021 Accepted: 16 July 2021 Published: 20 July 2021

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

as polyhydroxyalkanoates (PHAs) and polyhydroxybutyrate (PHB) or polylactic acid (PLA) [5–7].

In 2020, bioplastics accounted for only 1.2 million tonnes, less than 1% of all plastics produced [3]. However, the market for these polymers has been growing for years and the trend is expected to increase, with production estimated to increase by 11.28% from 2019 to the end of 2025 [8]. Specifically, PLA is one of the biopolymers with the greatest potential for industrial use as a substitute for petroleum-based polymers [9,10]. This is due to its great potential due to the combination of its mechanical properties, easy processability and its price compared to other biopolymers [11,12]. PLA has a suitable thermal stability and resistance to being industrially processed by injection molding, welding, thermoforming or extrusion [13]. Moreover, as it is biodegradable, one of its most interesting applications is mainly in the packaging sector [14]. In terms of mechanical properties, this material is comparable to non-degradable polymers in the "commodity" range [14,15]. However, there is a major limitation of PLA, its brittleness, which is a major drawback in the packaging industry [16]. To alleviate this disadvantage, there are several proposals. One solution to improve the flexibility of this polymer is its blending with other polymers such as polyethylene glycol (PEG), thermoplastic starch (TPS) or polybutylene succinate-co-adipate (PBSA), among others. However, the lack of miscibility between the components makes it difficult to achieve this improvement in toughness [5,11,17,18]. For this reason, some authors proposed modified vegetable oils (MVOs) as compatibilizing agents or even as plasticizing additives [12,19]. The reason for this is that they are a renewable and sustainable substitute for synthetic modifiers [20,21], which are also respectful of human health due to their non-toxicity, as they do not generate the migration of substances such as Bisphenol A (BPA), as is the case with conventional epoxy resins used in consumer products [22].

Several techniques can be used for chemical modification, such as epoxidation [23,24], maleinization [5,10,25], acrylation or hydroxylation. The MVOs available on the market today are mainly epoxidized soybean oil (ESBO) [25] and epoxidized linseed oil (ELO) [26]. Authors such as Garcia-Garcia et al. or Chieng et al. [26,27] have reported the improvement of PLA stiffness properties with these oils. Additionally, other authors investigated other modified oils for PLA formulations such as epoxidized cottonseed oil (ECSO) [24] or epoxidized palm oil (EPO) [27]. Another commercial option is maleinized linseed oil, which has provided excellent properties to PLA, as reported Ferri et al. [28].

In this work, the chemical modification of process of oils carried out is the maleinization. Maleinization is a chemical process, usually carried out in a single step, which consists of incorporating maleic anhydride molecules into the triglycerides that make up vegetable oils when conjugated carbon–carbon double bonds are present. For this purpose, there are several methods that can be used, such as the so-called "ene" reaction, Diels-Alder addition and free radical copolymerization, the first being the most favorable and the one used in this work [5,29]. This procedure requires a temperature of about 200 ◦C to result in the addition of the anhydric groups at the allylic positions of the fatty acid, as shown in Figure 1.

Both Brazil nut (*Bertholletia excelsa*) and hemp seed (*Cannabis sativa* L.) are interesting as MVOs, since both have an interesting lipid profile for functionalization. The Brazil nut is a brown fruit generally cultivated in the Amazon [30]. Its oil is an interesting object of study, since it has between 60.8% and 72.5% lipids [30,31]. Among them, it has 75.6% unsaturated fatty acids (UFA) [30], and both monounsaturated (MUFA) and polyunsaturated (PUFA), which means a significant amount of double bonds that allow its functionalization by chemical processes such as maleinization. On the other hand, hemp seed oil is also attractive for this purpose, as it is high in linoleic acid (55.3%) and linolenic acid (20.3%) [32]. Both allow their oil to be extracted by cold pressing with a good yield and also have good oxidative stability [31,33]. Thereupon, this work evaluates the potential of MBNO and MHO as a bio-based plasticizer to improve the ductile properties of PLA

related to its brittleness and the comparison of these with commercial maleinized linseed oil (MLO). MBNO and MHO as a bio-based plasticizer to improve the ductile properties of PLA related to its brittleness and the comparison of these with commercial maleinized linseed oil (MLO).

also have good oxidative stability [31,33]. Thereupon, this work evaluates the potential of

*Polymers* **2021**, *13*, x FOR PEER REVIEW 3 of 18

**Figure 1.** Schematic representation of the maleinization process of the triglyceride presents in vegetable oil. **Figure 1.** Schematic representation of the maleinization process of the triglyceride presents in vegetable oil.

#### **2. Materials and Methods 2. Materials and Methods**

LLC (Minnetonka, MN, USA) in pellet form.


*2.3. Bionic Microfluidic Device Design and Visualization of Injected Erythrocyte Trajectories*

The bionic microfluidic channel was designed using AutoCAD software and fabricated through standard lithographic method [17]. The microfluidic device was composed of polydimethylsiloxane (PDMS) and a glass plate, which is one layer with the height of 20 μm, and the detailed channel parameters are shown in Figure 1. After exposed to stable shear stress on a rheometer platform for 10 min, the whole blood sample was collected and injected into the inlet hole of the microchannel with the flow rate of 0.12 mL/min. The visualization of injected flowing erythrocyte trajectories in the microchannel were captured using a high-speed camera (M710, Phantom Co., New Jersey, USA) under the bright field provided by an inverted microscope (Leica DME6000 B, Leica Co., Wetzlar, Germany). The video analysis was processed using ImageJ software.


**5. Epilepsy of infancy with migrating focal seizures (EIMFS)**

### **5.1 Overview**

EIMFS, previously called malignant migrating partial seizures of infancy, is a rare and severe condition describe in 1995 [63]. EIMFS is characterized by focal "migrating" or "random" seizures beginning within the first 6 months of life, severe global developmental delay, and acquired microcephaly. Epilepsy is highly pharmacoresistant. At onset brain MRI is typically normal, but later it may show delayed myelination, thin corpus callosum, and cerebral atrophy [64]. Patients with intractable seizures have a progressive deterioration with major axial and limb hypotonia, loss of visual contact, and loss of other motor and social skills. Pyramidal and/or extrapyramidal features with athetotic movements may appear; about 18% of the patients died [65].

De novo KCNT1 mutations have been reported in about 50% of patients with sporadic EIMFS [64–66].

### **5.2 Seizures: symptoms and semiology**

It is accepted that the natural history of the EIMFS goes through three distinct phases [63]. The first phase starts in the first 6 months of life and lasts a few weeks or months; patients have sporadic seizures with frequency every few weeks or months. Seizures used to be focal motor with quick generalization or associated with autonomic features like apnea, flushing, or cyanosis [63, 67]. The second phase, called "stormy phase," arises between 1 and 12 months; seizures become more frequent, occurring in clusters several times a day or being practically continuous for several days. Seizure consists of lateral deviation of the head and eyes, twitches of the eyelids, unilateral clonic or tonic jerks of one or both limbs, apnea, flushing and/or cyanosis of the face, chewing movements, and secondary tonicclonic seizures [68]. Additionally, clinical manifestations may be "subtle" or absent even with the long duration of seizures.

The age at onset of the third phase is variable, from the end of the first year to the fifth year of age. This phase is typically seizure-free, even if interposing illnesses can trigger recurrent seizures or SS. Some patients can evolve to a WS [63, 65].

Migrating focal seizures (MFS) are seizures which occur usually in clusters that last for a few days and are then followed by a few weeks or months of recovery. Within a cluster, seizures are very frequent and may even extend to SS. Clusters increase in frequency within the first 2 years of life.

**89**

*Epileptic Encephalopathies in Infants and Children DOI: http://dx.doi.org/10.5772/intechopen.85378*

EEG is usually normal during the first months apart from slowing for many hours after long-lasting seizures. As the disease progresses, background activities become gradually diffusely slower with decline of physiological features. Activity may show alternating asymmetries, with slow activity that change from one hemisphere to another [63]. In seizure-free periods, sleep and wakefulness are obviously differentiated; nevertheless sleep spindles are rare, asynchronous,

Interictal abnormalities are usually absent at onset; spikes rapidly grow in frequency and develop multifocal within a few months; multifocal spike-and-wave

Epileptic discharges (ED) sequentially involve different areas of the brain, such as describing a random migration, without a specific pattern. The ictal pattern is characterized by rhythmic monomorphic activity in the alpha-theta frequency range, although delta waves, spikes, and spike-waves can also be observed. It is common for epileptic activity to remain limited to one region for a period of time and then decrease in frequency until stop, with a tendency to progressively involve an adjacent area. ED is continued with slow postictal activity without prolonged voltage decrement [63, 68]. When epileptic seizures are frequent, ED changes from one region to another and from one hemisphere to another so that consecutive focal epileptic discharges overlap resulting in a continuous and changing multifocal ictal activity and a very complex epileptic status pattern [11].

Migrating partial seizures are usually refractory to pharmacologic treatment though some cases have responded to bromide (60–80 mg/kg/day), with a termination of the seizures for several months after almost 3 weeks of therapy [69]. Even so, one should be aware of a potential bromoderma tuberosum, which could be appearing with high doses of potassium bromide therapy [70]. Intravenous levetiracetam (60 mg/kg) rapidly interrupted migrating partial status in two children with a good tolerability and safety [71]. Other successful treatments include a combination of sodium bromide, stiripentol, and levetiracetam [72], rufinamide and acetazolamide

**6. Epilepsy with myoclonic-atonic seizures (EMAS) or Doose syndrome**

In 1970, Herman Doose reported seizures in 51 previously normal children between 1 and 5 years of age described as myoclonic and astatic, frequently combined with absences and GTCS and tonic seizures [75]. Doose suggested a genetic etiology [76] and later refined his criteria and emphasized that tonic seizures are rare [77].

[73, 74], and stiripentol associated with clonazepam.

activities do not show any specific pattern and are not activated in sleep.

**5.3 Electroencephalography**

*5.3.1 Background*

and asymmetric.

*5.3.3 Ictal EEG*

**5.4 Treatment**

**6.1 Overview**

*5.3.2 Interictal abnormalities*