#### Abstract

The distribution characteristics of hydraulic gradient in embankment are closely related to seepage failure. Seepage failures such as flowing soil and piping will lead to serious damage and even the overall failure of embankment. The hydraulic conductivity has strong spatial variability, which changes the distribution of hydraulic gradient in embankment and increases the difficulty for predicting the embankment seepage instability. In this study, the distribution of soil hydraulic conductivity in a section of Shijiu Lake embankment was obtained by the permeability test. Based on Local Average Subdivision technique, a three-dimensional multilayer random field of embankment hydraulic conductivity was generated. Then, the mean and standard deviation of overflow point height and hydraulic gradient were calculated by the Monte Carlo method, which combined the generated three-dimensional random model and the deterministic analysis method of seepage field. Finally, the coefficient of variation (COV) of hydraulic conductivity (0.1, 0.3, 0.5, 0.7, 1.0, 2.0, and 3.0), the fluctuation scale in vertical direction (3 m) and the fluctuation scale in horizontal plane (3 m, 6 m, 12 m, 24 m, 36 m, and 48 m) were selected respectively for analyzing the random characteristics of embankment overflow point height and hydraulic gradient under the influence of different COV and fluctuation scale of embankment soil hydraulic conductivity.

#### 1. Introduction

There is strong randomness in the hydraulic conductivity of embankment soil, resulting in uncertainty in the calculation of seepage field. Different from the result of deterministic analysis, the responses of free surface, head and hydraulic gradient of seepage field change significantly, which is more consistent with engineering practice. Uncertainty analysis is closely related to random field theory, common random field generation methods include Fast Fourier Transform Method (FFT) [1], Turning Bands Method (TBM) [2, 3] and Local Average Subdivision Method (LAS) [4, 5]. The FFT method has the virtue of the high computational efficiency and calculates the effective random field quickly, while there is symmetry in the covariance in discrete random field, which easily causes significant error in some cases. TBM method discretizes the random field by setting the number and orientation of basic lines. Large quantity of lines results in high covariance accuracy, but huge calculation, which is difficult to meet the requirements of engineering application. The LAS technique for generating random field has the characteristics of high calculation accuracy, fast convergence, stable covariance calculation accuracy and small error in all cases, especially for strong randomness.

The spatial distribution of total water head, hydraulic gradient and flow rate are a typical feature of random seepage field, which is prerequisite for reliability analysis. Stochastic finite-element method is a common uncertainty analysis method. With its development, many different theories are proposed, and their applicable conditions and scope are also different. Taylor expansion stochastic finite-element method [6, 7] uses Taylor series expansion for the random variables at the mean point. The order of Taylor expansion is positively related to calculation accuracy but negatively related to the calculation efficiency. The first-order Taylor expansion method is commonly used. Perturbation expansion stochastic finite-element method [8, 9] performs first-order or second-order perturbation expansion on random variables. The calculation efficiency and accuracy are closely related to the number of perturbation expansion. Newman expansion stochastic finite-element method [10] expands random variables based on Newman series. The calculation efficiency of this method is relatively high and the accuracy meets the engineering requirements. However, the above methods ignore the high-order term and error term of the expansion polynomial, which limits the selection range of the COV of random variables and not suitable for the case with strong randomness. The increase of the order expansion can solve this problem, but significantly reduces calculation efficiency. Obviously, these characteristics restrict the application of these methods in engineering. Dam seepage is free surface seepage. Freeze [11] conducted stability and instability analysis on Earth dam and calculated the distribution of free surface. Bathe and Khoshgoftaar [12] gave the stochastic finite-element solution to the free surface seepage problem. Fenton and Griffiths [13], Ahmed [14, 15] believed that the hydraulic conductivity obeys lognormal distribution, and analyzed the free surface seepage of gravity dam. The results show that the height of overflow point by random analysis is limited compared with that in deterministic analysis. Zhou et al. [16] studied the seepage problem in heterogeneous porous media by the high-order numerical manifold method and got similar conclusions. Soil cracks introduce new seepage channels, which changes the distribution of hydraulic gradient [17–21]. Hydraulic gradient changes the stress and displacement field of soil, and has an impact on the stability of the structure [22–25].

These studies hardly consider that the soil of embankment is three-dimensional and multilayered, and the conclusions are far from the engineering practice. In this paper, considering the strong randomness of embankment soil hydraulic conductivity, a three-dimensional multimedia permeability random field of embankment based on LAS technique is generated, then, Monte Carlo simulation method is used to study the random characteristics of overflow point height and hydraulic gradient in embankment.

#### 2. Random Model of Hydraulic Conductivity

Seepage instability of embankment is a process of interaction between hydraulic gradient, soil strength and impermeability. As a complex multilayer material, the hydraulic conductivity of embankment soil has significant spatial variability. This strong spatial randomness has a significant impact on the random distribution of hydraulic gradient, and then affects the seepage stability of embankment. In order to obtain the random distribution of hydraulic conductivity of embankment soil under spatio-temporal evolution, many scholars have carried out a series of studies. Freeze [26] regarded clay as a homogeneous isotropic medium, and its hydraulic conductivity follows normal distribution or lognormal distribution in space. Hoeksema and Kitanidis [27], Sudicky [28] and Huang et al. [29] considered that the probability density function of soil hydraulic conductivity obeys lognormal distribution in space. In these theories, it is assumed that the hydraulic conductivity follows a certain mathematical distribution in space, and the regression analysis is carried out, then the spatial random characteristics of hydraulic conductivity is obtained.

Assuming that the hydraulic conductivity follows a lognormal distribution in space as shown in equation (1), the hydraulic conductivity of saturated seepage element can be obtained bywhere is the mean value of the hydraulic conductivity of a random field, is the variance of the hydraulic conductivity, is the mean value of logarithmic hydraulic conductivity , is the variance of . is the hydraulic conductivity assigned to the *i*th element in the random field, is the local average of a standard Gaussian random field, , over the domain of the *i*th element.

The random field correlation function is exponential, as shown inwhere *R* is the fluctuation scale of hydraulic conductivity in a certain direction and *r* is the distance between two points in the random field, always taking a positive value. Generally speaking, the scale of fluctuation is the distance over which points in the field are significantly correlated. It can be seen from equation (3) that when the fluctuation scale of the random field *R* increases, the value of the correlation function increases, which means that the correlation between points in the random field increases and the influence is greater. When the distance *r* between two points increases in the random field, the value of correlation function decreases and the influence between points will gradually weaken.

#### 3. Numerical Calculation Model

It is assumed that the embankment model consists of three parts, from top to bottom are embankment, layer 1 and layer 2. The dimensions of each part are shown in Figure 1, the embankment is 8 m high, the top of the embankment is 6 m wide and the bottom of the embankment is 54 m wide. The toe of the upstream slope is 30 m away from the left boundary of the model, and the toe of downstream slope is 60 m away from the right boundary of the model. Layer 1 is 4 m high in *Z* direction and 144 m long in *X* direction, while layer 2 is 20 m high in *Z* direction and 144 m long in *X* direction. When dividing the model into elements, considering the influence of fluctuation scale, the element size in the embankment is determined to be 3 m × 3 m × 1 m in the *X*, *Y*, *Z* direction respectively, layer 1 is determined to be 3 m × 3 m × 1 m in the *X*, *Y*, *Z* direction and layer 2 is determined to be 3 m × 3 m × 3.33 m in the *X*, *Y*, *Z* direction respectively. The total number of elements is 7488 and the total number of nodes is 8983 in the model.

The COV and fluctuation scale are closely related to the random field structure. The COV reflects the degree of random variables deviating from their mean value, the scale of fluctuation is the distance over which points in the field are significantly correlated. When generating the random field, it is assumed that the hydraulic conductivity is independent in the horizontal and vertical directions, firstly, the horizontal hydraulic conductivity is discretized randomly in a cuboid space, and then the vertical hydraulic conductivity is discretized randomly. After that, the anisotropic hydraulic conductivity of each element in the random field is obtained. According to the centroid position of the element, the hydraulic conductivity of the random field is mapped into the embankment model, and the hydraulic conductivity of each element in the embankment model can be obtained. Table 1 lists the mean of hydraulic conductivity at different parts of the model.

Figure 2 shows the three-dimensional multilayer random field of hydraulic conductivity in the model, with the COV of hydraulic conductivity is 0.3, the vertical fluctuation scale is 3 m, and the horizontal fluctuation scale is 6 m. Different gray levels are used to represent the hydraulic conductivity of the element, the darker the color, the smaller the hydraulic conductivity, on the contrary, the lighter the color, the greater the hydraulic conductivity. It can be seen that the gray levels of different elements vary greatly, When the water flows through the element with small hydraulic conductivity, the potential energy decreases rapidly, in contrast, when it flows through the element with large hydraulic conductivity, the water head decreases slightly.

In this study, the hydraulic conductivity is regarded as a random field variable which obeying lognormal distribution, and the COV is taken as 0.1, 0.3, 0.5, 0.7, 1.0, 2.0 and 3.0 respectively. The vertical fluctuation scale is fixed to 3 m, and the horizontal fluctuation scale is 3 m, 6 m, 12 m, 24 m, 36 m and 48 m. The Anisotropy ratio *ζ* can be obtained usingwhere is the scale of fluctuation in the horizontal plane and is the scale of fluctuation in the vertical direction.

The boundary conditions are as follows: the upstream water level is 31 m, and the downstream water level is 24 m, the bottom of the model and the boundaries on both sides are impervious. A realization refers to a single generation of the hydraulic conductivity random field and the subsequent finite-element analysis of the seepage field with ABAQUS. A Monte Carlo process involves a large number of realizations that eventually enable statistical statements to be made about the output quantities of the three-dimensional seepage. After 2000 independent realizations, the results are statistically analyzed, and the variation law of three-dimensional random seepage field of embankment is obtained.

#### 4. Results

##### 4.1. Overflow Point Height

Figure 3 shows the overflow point height calculated by ABAQUS in a single realization when the COV was 0.5, the vertical fluctuation scale was 3 m, and the horizontal fluctuation scale was 6 m. The water level of the upstream slope is in the same plane, and the height of the overflow point in the downstream slope shows uneven fluctuation. In the region with high overflow point height, the slope of the free surface is gentle, and the hydraulic gradient in the seepage path is relatively small, which is not easy to damage. Correspondingly, when the overflow point is low, the hydraulic gradient in the seepage path is high, and observation is required to ensure safety. The spatial variability of hydraulic conductivity leads to different overflow point height and hydraulic gradients in different sections of the embankment. Compared with deterministic analysis, it has advantages in describing the potential risk of embankment, and has better guiding significance for engineering practice.

Figure 4 shows the contour map of mean of the hydraulic gradient in the model, in the region far away from the upstream and downstream slope, the contour map of hydraulic gradient is nearly parallel. The reason is that in the upstream region, the boundary conditions are nearly the same everywhere, and the main seepage direction is approximately vertical downward, and in the downstream region, the main seepage direction is approximately vertical upward. The maximum hydraulic gradient appears below the overflow point in the downstream slope and close to the embankment toe, which is similar to the results obtained by the deterministic method. Previous studies show that this is also a dangerous part in theory and practice. At the same time, it verifies the correctness of the calculation method used in this study [30, 31].

Figure 5 shows the contour map of standard deviation of the hydraulic gradient in the model. Near the slope toe on both sides of the embankment and in layer 1, the contour map of standard deviation of hydraulic gradient is densely distributed, as well as, the corresponding contour map of hydraulic gradient is also densely distributed. This phenomenon is mainly due to the fact that the hydraulic conductivity of soil layer 1 is much lower than that of soil layer 2. The total head of seepage field decreases much faster in the soil layer with small hydraulic conductivity than in the soil layer with high hydraulic conductivity, At the same time, the hydraulic gradient generally reflects the change of water head per element length, Therefore, the hydraulic gradient of layer 1 is high, while that of layer 2 with high hydraulic conductivity is smaller, this phenomenon accords with the general law of seepage field.

Figure 6 shows the variation of overflow point height with the hydraulic conductivity for different values of COV. The height of overflow point in downstream slope decreases gradually during the gradual increase of COV from 0.1 to 3.0. When , the slope of the curve is large, and the height of the overflow point decreases with a faster rate, When , the slope of the curve decreases gradually, and the reduction rate of the overflow point height slows down, When , the reduction rate decreases significantly. Therefore, it can be concluded that when the mean of hydraulic conductivity is a constant, the overflow point height is inversely proportional to the COV of the random field variable. The main reason is that the spatial variability of hydraulic conductivity leads to the seepage path tortuous and the streamline longer.

When the COV begins to increase, the variation of the hydraulic conductivity in the random field increases, and its value range increases accordingly. At this time, the probability of obtaining a smaller or larger value of hydraulic conductivity increases. According to the principle of minimum potential energy, when flowing through the region with large difference in hydraulic conductivity, the occurrence of seepage always tends to avoid the region with small hydraulic conductivity and flow through the region with large hydraulic conductivity, at this time, the streamline begins to appear twists and turns. Considering that the seepage always flows from the high head to the low head, the height of the overflow point is reduced. When the COV increases further, the probability of obtaining a smaller or larger value of the hydraulic conductivity increases accordingly, resulting in a further reduction of the overflow point height. When the seepage occurs in the region with small hydraulic conductivity, the energy loss increases and the water head decreases greatly when the seepage flows through the unit length. The equipotential surface decreases greatly compared with the uniform medium, resulting in the reduction of the overflow point height.

Figure 7 shows the variation of the overflow point height with the anisotropy ratio for different values. It can be seen that the mean of overflow point height shows a slight growth trend with the increase of fluctuation scale. When a deterministic analysis is adopted, the overflow point height is a constant, which means the maximum value of the height. When the coefficient of variation approaches 0, the soil can be regarded as homogeneous, and the height curve is closer to a straight line. When the coefficient of variation increases, the hydraulic conductivity is different from that of homogeneous soil obviously, the height of overflow point decreases, and the curve is no longer an approximate straight line. When the COV of hydraulic conductivity is small, the curve approaches a straight line. When the COV is greater than 0.7, the curve shows an upward trend obviously. The main reason is that when the COV is small, the soil is relatively uniform and the influence of fluctuation scale is weak, with the increase of COV, the soil heterogeneity gradually increases, a large fluctuation scale means that the hydraulic conductivity is interrelated in a large range, and resulting in a certain consistency within this range, which reduces the soil heterogeneity. At this time, the fluctuation scale has an obvious impact on the height of the overflow point, resulting in an upward trend of the curve.

when , the minimum value of the curve appears while the anisotropy ratio is equal to 2, which shows that within a certain range of COV, the minimum height of overflow point does not appear at the minimum anisotropy ratio, but somewhere near the minimum. This phenomenon plays an important role in statistical the height of overflow points, which may cause other errors in the project further and have a negative impact on the feasibility, design and construction of the project.

##### 4.2. Hydraulic Gradient

Figure 8 shows the variation of the mean of hydraulic gradient with the hydraulic conductivity for different values of COV at the toe of the downstream slope in the model. In general, the hydraulic gradient is directly proportional to the COV, When , the mean of hydraulic gradient increases slowly with the increase of COV, When , the mean of hydraulic gradient increases more rapidly with the increase of COV. The reason for this phenomenon is that when the COV is small, the soil in embankment is relatively uniform, and there is little difference between the random solution and the deterministic solution of the hydraulic gradient in embankment. With the increase of COV, the difference of hydraulic conductivity of each element increases, in the soil element with high hydraulic conductivity, the hydraulic gradient is obviously small due to the little water retaining effect of soil, while in the element with small hydraulic conductivity, the hydraulic gradient increases obviously.

Figure 9 shows the variation of the mean of hydraulic gradient with the anisotropy ratio for different values at the toe of the downstream slope, When , the curve drops sharply, and when , the curve is relatively flat. This shows that when the COV is large, the smaller anisotropy ratio has a great impact on the hydraulic gradient, and the mean of the hydraulic gradient decreases rapidly. With the increase of anisotropy ratio, the reduction ratio of hydraulic gradient becomes slow. The author believes that when the anisotropy ratio increases, a more uniform horizontal seepage channel is formed locally in the soil, water tends to pass through the horizontal seepage channel rather than downward seepage, resulting in the rapid reduction of the hydraulic gradient. With the further increase of anisotropy ratio, a complete horizontal seepage channel is formed, and the change of hydraulic gradient tends to be slow [32].

Figure 10 shows the variation of the standard deviation of hydraulic gradient with the hydraulic conductivity for different values of COV. With the gradual increase of the COV, the standard deviation of hydraulic gradient shows a rapid increasing trend, and when the fluctuation scale is small, the curve is located in the upper part in Figure 10, and the location of the curve decreases step by step with the increase of the fluctuation scale. The reason is that when the fluctuation scale is large, the soil is in a strong correlation state in the model, and the hydraulic conductivity of soil tends to be uniform in the model, so as to reduce the standard deviation of hydraulic gradient.

Figure 11 shows the variation of the standard deviation of hydraulic gradient with the anisotropy ratio for different values. The standard deviation of hydraulic gradient decreases slowly with the increase of anisotropy ratio. When the COV is small, the standard deviation of hydraulic gradient hardly changes, indicating that the standard deviation is not sensitive to the change of anisotropy ratio. Therefore, in practical engineering or searching for dangerous sections in embankment, the influence of COV on the standard deviation of hydraulic conductivity should be considered, and the anisotropy ratio should be considered when necessary.

#### 5. Conclusion

This study focuses on the layered structure of embankment soil and the strong variability of hydraulic conductivity. A three-dimensional multilayer random field of hydraulic conductivity is generated, and a Monte Carlo process involves 2000 realizations are used to solve the overflow point height and hydraulic gradient in embankment. Based on this study, the following conclusions can be drawn:(1)The height of the overflow point in the downstream slope shows uneven fluctuation. When the overflow point is high, the hydraulic gradient in the seepage path is small and has no damage effect. In contrast, the hydraulic gradient on the seepage path is large with low overflow point, which needs to be observed.(2)The maximum hydraulic gradient appears below the overflow point of downstream slope and near the embankment toe, where the maximum standard deviation of hydraulic gradient also appears. It means that the hydraulic gradient at these positions is strong, highly volatile, and easy to be damaged.(3)When the COV of hydraulic conductivity increases, the height of the overflow point decreases gradually, while the height of the overflow point increases slowly with the increase of the anisotropy ratio.(4)When the COV increases, the hydraulic gradient goes up obviously, which effects embankment instability. the anisotropy ratio increases with the hydraulic gradient decreasing significantly, which reduces embankment instability.(5)When the COV increases, the standard deviation of hydraulic gradient inclines significantly, which increases the fluctuation of hydraulic gradient and expands the unstable factors. the increase of anisotropy ratio helps the standard deviation decrease slowly, which also reduces the risk to a certain extent.

#### Data Availability

The datasets generated during the current study are available from the corresponding author on reasonable request.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Authors’ Contributions

The manuscript was approved by all authors for publication.

#### Acknowledgments

The authors would like to acknowledge financial support from Key Scientific Research Projects of Colleges and Universities in Henan Province under Grant number 22B570002 and Henan College Students’ Innovation and Entrepreneurship Training Program under Grant number S202110480034.