Numerical simulation of the free surface and water inflow of a slope, considering the nonlinear flow properties of gravel layers: a case study

Groundwater is an important factor of slope stability, and 90% of slope failures are related to the influence of groundwater. In the past, free surface calculations and the prediction of water inflow were based on Darcy's law. However, Darcy's law for steady fluid flow is a special case of non-Darcy flow, and many types of non-Darcy flows occur in practical engineering applications. In this paper, based on the experimental results of laboratory water seepage tests, the seepage state of each soil layer in the open-pit slope of the Yanshan Iron Mine, China, were determined, and the seepage parameters were obtained. The seepage behaviour in the silt layer, fine sand layer, silty clay layer and gravelly clay layer followed the traditional Darcy law, while the gravel layers showed clear nonlinear characteristics. The permeability increases exponentially and the non-Darcy coefficient decreases exponentially with an increase in porosity, and the relation among the permeability, the porosity and the non-Darcy coefficient is investigated. A coupled mathematical model is established for two flow fields, on the basis of Darcy flow in the low-permeability layers and Forchheimer flow in the high-permeability layers. In addition, the effect of the seepage in the slope on the transition from Darcy flow to Forchheimer flow was considered. Then, a numerical simulation was conducted by using finite-element software (FELAC 2.2). The results indicate that the free surface calculated by the Darcy–Forchheimer model is in good agreement with the in situ measurements; however, there is an evident deviation of the simulation results from the measured data when the Darcy model is used. Through a parameter sensitivity analysis of the gravel layers, it can be found that the height of the overflow point and the water inflow calculated by the Darcy–Forchheimer model are consistently less than those of the Darcy model, and the discrepancy between these two models increases as the permeability increases. The necessity of adopting the Darcy–Forchheimer model was explained. The Darcy–Forchheimer model would be applicable in slope engineering applications with highly permeable rock.


Introduction
In China, from 1949 to 2011, even the most conservative estimates of the deaths from landslides exceed 25 000 people, which was a death toll of more than 400 people each year, and the average annual pecuniary loss was $50 million dollars [1]. Groundwater is one of the main factors causing instability and sliding failure of the soil slope. A total of 90% of slope failures, for either natural slopes or artificial slopes, are related to groundwater [2]: more specifically, porewater pressures and seepage forces. It can be seen that understanding the seepage mechanism of slopes is necessary to prevent the occurrence of landslides.
The effect of water on the slope stability can be seen by considering the slope with different positions of the free surface [3]. Slope failures associated with the water level fluctuations were reported by many researchers, such as Lane [3], Jappelli & Musso [4], Shi & Zheng. [5], Jia et al. [6], Wang et al. [7] and Lu et al. [8]. The position of the free surface has a great effect on the slope stability; the stability safety coefficient of a soil slope gradually decreases with the continuous increase in the free surface area, until slope sliding occurs [9]. When the velocity of flowing reservoir water is uniform, a simplified formula, based on Boussinesq's theory, can be used to calculate the location of the free surface; a numerical analysis indicated that when the free surface is one-third of the slip body, the stability coefficient of the slope is minimal [5]. To investigate the influence of the drawdown on the slope stability and pore water pressure, Wang et al. [7] performed some laboratory physical model tests. Lu et al. [8] adopted an appropriate method to evaluate slope stability during earthquakes, including the effect of the groundwater level; this method extended the scheme of the limit equilibrium method. The free surface fluctuations may induce the change in the pore water pressure of the soil slope. Therefore, calculating the free surface of the soil slope is very useful to analyse the stability of the slope.
Moreover, it is well known that the water inflow of a slope has a very strong influence on the slope stability. Some fine particles are washed away due to the effect of water seepage, which is a main cause of mass loss. Hu et al. [10] adopted a group of orthogonal tests to investigate the influence of water capacity on sand-soil slope stability. To study the water inflow of an open-pit slope of the Ekou Iron Mine in China, Liu et al. [11] built a numerical seepage model based on GeoStudio and AutoCAD-3Dmine. Slope stability analyses and supporting studies under the condition of a larger foundation pit water inflow have been studied by Ma [12].
The above research all assumes that the water seepage in the slope obeys Darcy's law. However, the slope of the Yanshan Iron Mine contains two high-porosity gravel layers, which are very sensitive to the flow of groundwater. The results of laboratory experiments indicate that the fluid migration in these gravel layers follows non-Darcy flow and not Darcy's law. Li [13] also performed similar studies.
Until now, many studies on the Darcy fluid flow in this slope have been carried out, but the analysis of the non-Darcy flow of the slope is less reported in relevant publications. According to different theories and methods, including Darcy theory and non-Darcy theory, the results of these studies are sometimes very different. To investigate this problem, this paper first expounds the general situation of the Yanshan Iron Mine in §2. In addition, this paper presents experimental equipment used to test the specimens of soil which form different layers of the Yanshan Iron Mine slope, and the seepage properties, including the non-Darcy coefficient, are determined from measurements and calculations, as described in §3. A coupled model that includes two flow fields is first established based on the theory of conservation of mass and pressure balance; this model couples the equations of the Darcy flow field and the Forchheimer field, as described in §4. Then, the numerical results, the effect of the porosity of the gravel layers on the water inflow and the free surface are described in §5. Finally, §6 provides conclusions and discusses the implication of these findings.    in the testing system, as shown in figure 4. In addition, figure 4d illustrates the system connections and principle: the water supply system (5) consists of a constant pressure water supply and constant flow water supply. The volume of the constant pressure pump is 12 l, and it can provide water pressure at 2 MPa to the upstream of the sample. The maximum flow during the tests is 25 l min −1 . In addition, the constant flow pump can provide a stable water flow at 0-3 l min −1 . The air compressor (6), pressure regulator (1), air cylinder (8) and plunger (9) are used to apply the axial loading pressure. The porous plates (11) can ensure that the fluid flows evenly and protects the structure of the porous media so that it is not destroyed under high water pressures. The water pressure on the upstream and downstream ends of the sample are measured by pressure sensor (7)    to the top, until the sample saturation was above 95%. In addition, the atmospheric valve (16) can be used to discharge air if there is air trapped in the tube. In the flow test, if the outlet (13) is directly open to the atmosphere, air fluctuation often causes fluid instability and the veracity of the measurement results cannot be ensured, even when using high-precision sensors. To avoid the outlet effects, back pressure must be applied at the outlet. The speed measuring device (15) includes an electronic balance (TCL 12000), in the range of 0 to 12 000 g ± 0.1% FS, and large beakers. Data collection and analysis are performed by the data acquisition system (3), and the sampling interval is 3 s.

Experimental materials
The soil specimens used in the tests were taken from the eastern slope of the Yanshan Iron Mine, China, and were prepared in the laboratory. The soil columns were sampled from a borehole.   boundary effect, before the tests, the inside of the test cylinder was coated with a thin layer of petroleum jelly. A vacuum pump was used to remove any air bubbles from the test cylinder. The measurements were carried out under a stable flow. The main method to experimentally measure the flow rate was the weighing method. Figure 5 shows the preparation process before the experiments were conducted.

Water flow experiments
To evaluate the permeability of each layer, all experimental data and the regression relationships between the hydraulic gradient and flow velocity were analysed and are depicted in figure 6. For the silt, fine sand, silty clay and gravelly clay layers, the relationships are linear, and the flow behaviour is Darcy flow. However, these relationships are nonlinear for the gravel layers; therefore, the flow behaviour is considered nonlinear flow [14]. The one-dimensional Darcy equation can be written as x (m -1 ) y (m -2 ) T1 y = 9.78 × 10 6 x + 4.69 × 10 9 R 2 = 0.99 T2 y = 7.18 × 10 6 x + 1.38 × 10 9 R 2 = 0.99 T3 y = 5.07 × 10 6 x + 5.72 × 10 8 R 2 = 0.99 T4 y = 3.81 × 10 6 x + 3.51 × 10 8 R 2 = 0.99 T5 y = 2.71 × 10 6 x + 3.07 × 10 8 R 2 = 0.99 T6 y = 1.91 × 10 6 x + 1.53 × 10 8 R 2 = 0.99 Figure 7. Relation between y and x. where J is the hydraulic gradient, μ is the dynamic viscosity, k is the permeability, g is the gravitational acceleration, ρ is the fluid density and v is the superficial velocity. The non-Darcy groundwater seepage problem is very complicated, and theoretical studies of this problem still use the Darcy flow theory. Until now, there was no formula that could describe this flow characteristic more accurately [15]. Through numerous tests, the empirical formulae of nonlinear flow were obtained. There are two basic forms: the Izbash equation and the Forchheimer equation [16]. The Forchheimer equation has a clear physical meaning and theoretical basis and is, therefore, widely used [17][18][19], the Forchheimer equation can be written as where ∇p is the pressure gradient and β is the non-Darcy coefficient.
To measure k and β, ∇p/μν is denoted as y and ρv/µ is denoted as x. Equation (3.2) then simplifies to According to the laboratory experimental data, in the x-y coordinate system, a straight line is defined with slope β and y-axis intercept 1/k. The test results of the soil permeability parameters are shown in table 1. Figure 7 shows the determination of k and β in the gravel layers with different porosities. The results of k and β and the porosity of the gravel layers are shown in table 2.
In practical engineering, porosity is easier than permeability to obtain experimentally. The Kozney-Carman (kc) equation is one of the most famous formulae for liquid or gas flow in porous media and has been used to calculate the permeability of porous media worldwide for a long time. Eventually, the permeability could be summarized as a function of the particle diameter and porosity [20]. Permeability can be defined as where d is the particle diameter and φ s is the shape coefficient.   In the gravel layers, the average size and shape coefficient of the particles are invariant; hence, d and φ s may be considered constants. Then, the relationship between permeability and porosity can be expressed as follows: where a is constant ( figure 8).
Through data fitting, we also present a new correlation: The correlation coefficient of equation (3.6) is 0.98; therefore, equation (3.6) is reasonable to use for practical applications in which the permeability of gravel layers needs to be determined. The porosity, permeability and tortuosity are probably the most important factors determining the value of the non-Darcy coefficient. However, determining a realistic and precise tortuosity is not easy. Therefore, many scholars use only the porosity, n, and permeability, k, to characterize the non-Darcy coefficient, β [21,22], which can be written as   The groundwater flow in the silt layer, fine sand layer, silty clay layer, gravelly clay layer and bedrock is mainly driven by the water pressure. The seepage velocity is very slow, the inertia force is considerably less than the viscous drag, and the flow state of the groundwater is laminar flow. The relation between the fluid velocity and pressure gradient are linear, consistent with the Darcy law at the macroscale, and can be expressed as where the subscript d represents the Darcy seepage field. Under unsteady-state conditions, the seepage continuity equation is where t is the time, f d is the source term, ∇ is the gradient operator and ∇· is the divergence operator.

Forchheimer flow in the gravel layers
In  (4.3) and the continuity equation is where the subscript f represents the Forchheimer seepage fields, c f is the acceleration coefficient and f f is the flow mass force. In addition, the permeability k f and the non-Darcy coefficient β f of porous materials considerably affect the water flow through the media. The relation between the permeability and non-Darcy coefficient is described by equation (3.9).

Boundary conditions for the transition between two flow regions
During the whole process of groundwater flow in the slope, both the flow velocity and the pressure are continuous. Therefore, the transition boundary conditions can be defined as The element birth and death method is used in the FELAC 2.2 software to solve the seepage flow problem with a free surface. A term is added to equation (4.9) so that when the water head is less than its position head, the left stiffness term becomes small by setting em = 1.0 −12 .

Numerical solution of the model
To use an iteration algorithm to calculate the free surface, setting an overflow point is necessary to balance the water head (h) with the position head (H). As the overflow point approaches, the boundary above the overflow point needs to be updated, and then the equations will be reorganized. Figure 10 illustrates the entire iterative procedure. In addition, according to the principle of virtual work, the Darcy equation and Forchheimer equation are represented in the form of a weak integral. To solve the model, the finite-element method is combined with a finite volume method. Then, the numerical programs are developed with FELAC 2.2 software. The specific process of solving the model has been illustrated in detail in the literature [23].

Numerical model
According to the geologic section of the Yanshan Iron Mine (figure 2), a two-dimensional numerical model is established, as shown in figure 11. The length of the numerical model is 140 m, and its height is 50 m. The computing area is divided into 37 076 four-node quadrilateral elements. According to the hydrological conditions, the boundary conditions of the model are as follows. The constant water heads, H = 45 m and H = 20 m, are set at the right boundary (red boundary), which is near the Xinhe River, and the left boundary (blue boundary) of the slope, respectively. At the left boundary, the other slope sides are directly exposed to the atmosphere; therefore, they are free boundaries. At the interface of the   The pressure loss is caused by the combination of the viscous resistance and inertial resistance, and more rapid decreases in pressure are observed when using the Forchheimer equation than when using the Darcy equation. Therefore, the results of the Darcy model are rather high.

Further discussion
As the water burst into the open pit, the permeability of the gravel layers increases with the loss of small particles due to the high-velocity water flow. Under certain water pressure conditions, the permeability and non-Darcy factor of the gravel layers not only directly affect the viscous drag and inertia resistance but also control the distribution of fluid pressure in the slope, the location of the free surface and the amount of water inflow. In addition, there is a one-to-one relationship between the permeability and non-Darcy factor [17,24]. Therefore, a sensitivity analysis of the permeability of the gravel layers is important to guide the accurate prediction of the amount of water inflow and the location of the free surface.

The effect of porosity on the overflow point
Determination of the seepage overflow point is very important in the analysis of soil-slope stabilization problems. To study the effect of the permeability of the gravel layers on the overflow point location, five cases are considered. Since some of the small particles are lost, the following five situations are considered: the porosity of the gravel layers are 0.1, 0.2, 0.3, 0.4 and 0.5, and the corresponding permeabilities are 2.28 × 10 −11 , 2.31 × 10 −10 , 1.02 × 10 −9 , 3.29 × 10 −9 and 9.25 × 10 −9 m 2 , respectively (calculated by using equation (3.6)). The non-Darcy factor, β, is calculated by equation (3.10). In addition, the permeability of the other layers remains unchanged. Table 4 shows the various seepage field physical quantities influenced by the values of n.
When the permeability of the gravel layers is low, the nonlinear phenomenon of the seepage flow is unclear, and the calculation results from both models are very similar. The height of the overflow point is not very different between the two models, therefore, the error of model 1 is relatively small. With an increase in the permeability, the height of the overflow points of model 1 and model 2 gradually   increase. Additionally, the discrepancy between the two results increases, and the height of the overflow point calculated by model 2 is consistently lower than that of model 1. When the permeability is higher, the water line is relatively gentle, while the relative curvature of the water level is greater due to the low permeability. In the seepage outlet area, the increase in the hydraulic gradient is obvious. Accordingly, it is necessary to calculate the free surface with model 2 to simulate large pores and high permeability. The phenomenon of non-Darcy seepage is less apparent when the flow is through less permeable porous media. However, the high permeability layers and low-permeability layers are mixed together in the natural geological environment or slope, complicating the seepage field. To more objectively and accurately determine the location of the overflow point, it is necessary to use a higher-precision calculation method, especially in the areas where the flow velocity is rapidly changing and in the areas that need to be studied in slope engineering.

The effect of porosity on the water inflow
The relation curves that describe the variation in water inflow with time are presented in figure 16.
The results show that the change in the amount of water inflow can be divided into three stages: (1) accelerated growth; (2) slow growth; and (3) steady flow. In the accelerated growth period, the rate of increase in the water inflow calculated by using the two models are fairly similar under the condition of low permeability, but they differ considerably with greater permeability, and the rates of the increase in the water inflow calculated by using model 1 is greater than that of model 2. The higher the permeability is, the greater the difference between the two model results. In addition, the amount of water inflow can become stable, even within a very short time, because in most cases, the transient fluid flow through the porous media decays very quickly, and the time derivative terms of equations (4.2)-(4.4) also decrease quickly, transforming the transient equations into steady-state equations [25]. According to our analysis, the results obtained through the two models were found to have few discrepancies in the steady flow period under the condition of low permeability. With the gradual increase in permeability, the maximum water inflow increases gradually and the maximum water inflow calculated by model 1 is consistently greater than that of model 2. That result is because the Forchheimer equation includes the square of the water velocity, considerably slower than the water velocity used in the Darcy equation under the same condition. Figure 17 shows that the maximum water inflow varies with permeability. Under certain boundary and initial conditions, the maximum water inflow increases with the increase in permeability, while the maximum water inflow of model 2 increases more slowly than that of model 1. Thus, it is necessary to adopt model 2 to predict the maximum water inflow for the open-pit slope that contains the highpermeability gravel layers. If we use model 1 to calculate the maximum water inflow, the results would greatly deviate from the actual situation. As the permeability increases, the difference between these two models increases.

Conclusion
The permeability of soil layers is obtained according to laboratory seepage experiments. The functions of the permeability versus porosity and non-Darcy factor versus porosity and the relationship among them were determined. A coupling model is established to study the water flow behaviour in the Yanshan Iron Mine. The effect of porosity on the overflow point and water inflow was discussed. The main conclusions are drawn as follows: 1. Laboratory tests show that the seepage behaviour in the silt layer, fine sand layer, silty clay layer and gravelly clay layer follow the traditional Darcy law, while the gravel layers show clear nonlinear characteristics. The permeability increases exponentially, and the non-Darcy coefficient decreases exponentially with an increase in the porosity of the gravel layers, and the relations among the permeability, the porosity and the non-Darcy coefficient are established as equations to use in future engineering applications. 2. For a slope containing gravel layers with high permeability, considering the seepage in the slope during the transition from Darcy flow to Forchheimer flow, a coupled mathematical model that considers both of these flow fields is established with Darcy flow in the low-permeability layers and Forchheimer flow in the high-permeability layers. This method provides a model base for predicting the water inflow and computing the free surface of a slope. 3. This model is used in the seepage analysis of a cross section of the eastern slope of the Yanshan Iron Mine. By comparing the two solutions from the Darcy model and the Darcy-Forchheimer model, the results indicate that the free surface calculated by the Darcy-Forchheimer model is in good agreement with the in situ measurements; however, there is an evident deviation of the simulation results from the measured data when the Darcy model is used. Therefore, the model is suitable for use in high-permeability slope engineering. 4. The sensitivity analysis of the gravel layers indicates that with the increase in the permeability of the gravel layers, the water inflow increases, and the location of the overflow point moves up the slope. The results of the Darcy model are greater than those of the Darcy-Forchheimer model, and the greater the permeability is, the greater the difference between the two models. Accordingly, it is necessary to calculate the free surface and water inflow with the Darcy-Forchheimer model when the slope includes large pores and high permeability. 5. In slope engineering, although the calculation results of the Darcy model are more conservative in terms of safety, the precision and objectivity of the description of the seepage behaviours can be improved by applying the Darcy-Forchheimer model. In mine drainage engineering, the more accurate the prediction of water inflow is, the lower the economic cost.