Study of Anti-Sliding Stability of a Dam Foundation Based on the Fracture Flow Method with 3D Discrete Element Code

Fractured seepage is an important factor affecting the interface stability of rock mass. It is closely related to fracture properties and hydraulic conditions. In this study, the law of seepage in a single fracture surface based on modified cubic law is described, and the three-dimensional discrete element method is used to simulate the dam foundation structure of the Capulin San Pablo (Costa Rica) hydropower station. The effect of construction joints and developed structure on dam stability is studied, and its permeability law and sliding stability are also evaluated. It is found that the hydraulic-mechanical coupling with strength reduction method in DEM is more appropriate to use to study the seepage-related problems of fractured rock mass, which considers practical conditions, such as the roughness of and the width of fracture. The strength reduction method provides a more accurate safety factor of dam when considering the deformation coordination with bedrocks. It is an important method with which to study the stability of seepage conditions in complex structures. The discrete method also provided an effective and reasonable way of determining seepage control measures.


Introduction
Seepage liquid is usually regarded as external load in a dam foundation or in rock engineering. The methods of evaluating the stability of complex infiltration load and of the effect of hydraulic-mechanical coupling are major challenges in civil engineering. A discrete fractured rock mass is rich in joint surfaces when compared with a porous medium. The directionality, connectivity, openness, filling, and roughness of joints contribute to the characteristics of inhomogeneity, anisotropy, and discontinuity in fractured rock mass, and significantly influence the seepage in a fractured rock mass. The above situation should be considered to determine whether the seepage characteristics of fractured rock are in accordance with the cubic law of seepage. Incorrect uses of the porous media seepage analysis method will easily lead to unreasonable results.
Traditionally, there are usually three types of seepage calculations in the analysis of fractured rock mass (Jing, 2003) [1]: the Continuum methods, Discontinuum methods, and a Hybrid continuum/discontinuum models. The equivalent porous media model as a continuum method is typically used and developed. He et al. (2011) [2] proposed a new approach for equivalent porous media model to evaluate equivalent permeability tensors, the finite element method was used in the study of single-hole packer tests. When it comes to the computational mechanics of discontinua and

Seepage in a Single Fracture
The flow of a fluid in a single fracture is usually described by the Navier-Stokes equation, which is based on momentum and mass conservation. The equation assumes that the fluid through the fracture is a stable laminar flow with a constant density and viscosity, and, moreover, that the fracture wall is impermeable. The vector form of the Navier-Stokes equation is represented as where ρ is the fluid density (kg/m 3 ), µ the fluid viscosity (N·s/m 2 ), → v = v x , v y , v z the fluid velocity vector (m/s), and p(x, y, z) the hydrodynamic pressure (Pa). However, the Navier-Stokes equations are difficult to solve using nonlinear partial differential equations. Therefore, a simplified treatment is usually used; namely, the local cubic law (LCL) that approximates a three-dimensional flow field into a two-dimensional flow field (zero velocity in the z direction) by the Stokes equation. The fracture surface is assumed to be rough, and the viscous force in the flow field plays an important role. The formula for the local cubic law is as follows: where b (m) is the local aperture parallel to the z axis, γ the bulk density (Pa), µ the fluid viscosity (N·s/m 2 ), and H (x, y) the average hydraulic head (m). The two-dimensional flow field of the LCL was approximated as a one-dimensional flow field in cubic law. The fracture surface is assumed to be a pair of parallel plates separated by a certain distance. Rock mass is assumed to be impermeable and seepage only occurs in the fracture surface. It is considered that the seepage rate of a single fracture is proportional to the third power of the fracture opening: where Q is the single wide flow rate (m 3 /s), γ the bulk density(Pa), W the crack width(m), µ the dynamic viscous coefficient (N·s/m 2 ), b the local aperture opening parallel to the z axis(m), and Engelder, Scholz, Raven, and some others proved through experiments that the cubic law is only suitable to describe a single fracture with a smooth wall, large opening, and no filling of the fracture. Tsang (1987) [18] pointed out that the infiltration in the fracture is actually a trench flow rather than a laminar flow because of the change in the fracture aperture size and the existence of nonconforming regions. Gudmundsson (2000) [19] proposed that the fracture aperture depends not only on the associated stress field, but also on its controlling dimension, thus realistic models applied to fluid flow in rock fractures should take it into consideration. In order to consider the effects of fracture roughness and aperture variation on permeability, Witherspoon (1983) [20] introduced the concept of the hydraulic conducting aperture and formed the generalized flow law as follows: where f is the joint roughness coefficient, C is a constant depending on the flow geometry and fluid properties. If n is taken as 3, Equation (5) degenerates into the cubic law. The hydraulic conducting aperture b h (µm) and the mechanical conducting aperture b (µm) have the following relationship: Barton and Bandis (1985) [21] proposed an empirical formula based on based on a constitutive model relating, mechanical aperture b E (µm), the theoretical smooth wall conducting aperture b e (µm), and joint roughness coefficient (JRC). The conductivity in units of length squared k H (m 2 ) can be estimated based on conversion from the b E to b e , which is nonlinear with b e .
Barton's model has a higher degree of agreement with the test results than some linear relationships (Rutqvist, 1996) [22]. Most experiments show that when the opening degree is reduced to a certain extent, the cubic law fails, and in this case the generalized cubic law (n = 3, f = 1) is still applicable. When the aperture size continues to decrease, the generalized cubic law is invalid. In this case, the generalized flow law should be adopted. Detournay (1980) [23] has proposed a relation between the hydraulic aperture, u h (µm), and the joint deformation, u m (µm): where f is a correction factor that reflects the influence of the roughness on the tortuosity of the flow, u m is defined as the average distance between the joint surfaces, and u m0 is the initial aperture. This relation permits a more general fitting of the "cubic law" to the experimental results. Assuming the following dependence of the joint conductivity (k H ) on hydraulic aperture, Alvarez et al. (1995) [24] reanalyzed some previously published data on joint conductivity, and have confirmed the validity of the "modified cubic law" (MCL) as follows: The MCL provides an approach for accurately calculating the hydraulic properties and local vertically integrated flow fields for rough and tortuous fractures. When it comes to the analysis of multifractured rock-mass seepage, mature three-dimensional, discrete element analysis software can deal with it well based on the MCL. On this basis, we discuss the anti-sliding stability of the dam foundation because it is suitable for integration with, and the improvement of, fracture network models.

Numerical Simulation Method and Verification
The presumptions made in simulating the seepage flow of a fractured rock mass are as follows: (a) The structural plane is permeable and the rock mass is impermeable; (b) structural plane seepage obeys the cubic law; and, (c) the hydraulic-mechanical coupling on the structural plane is considered, and the seepage rate of the structural plane under different stress conditions will change with the fracture width.
Cappa (2005) [25] conducted a series of in situ confirmatory tests of fracture seepage in the Southern French Alps. The site consists of impermeable marl at the bottom and upper 15 m accumulation of fractured limestone, as shown in Figure 1a. The spring water stored by the fault appears to be perpendicular to the nonpermeable fault. A 10 m thick layer of waterproof concrete was installed on the surface to avoid water leakage at the point of discontinuity, and thus the site can be regarded as a 30 × 30 m 2 natural reservoir. Water pressure can be controlled manually by a sluice, and is maintained at a level 10 m above the sluice when the sluice is closed; the stored water will flow out naturally when the sluice is opened.  Bedding plane S2 7.6 × 10 −6 3.4 × 10 −11 9.6 × 10 −5 20 -8 Bedding plane S2 3..8 × 10 −6 2.1 × 10 −11 6.8 × 10 −5 20 -9a Bedding plane S2 3.2 × 10 −6 1.8 × 10 −11 6.2 × 10 −5 20 -9b Bedding plane S3 9.0 × 10 −7 0.5 × 10 −11 3.3 × 10 −5 20 -10 Fault Piezometers were arranged densely in the test area, and six sets of stress-strain measurements were acquired at three discontinuous measuring points. The water injection test of a fault (designated F12) was carried out and the variation of its width during the water-filling process was measured. The hydraulic and mechanical properties of the main fault structures are shown in Table 1. The hydraulic-mechanical coupling of a fractured rock mass is complex, and it is difficult to obtain accurate solutions for it. Moreover, the development of computer numerical methods provides an important method of hydraulic coupling analysis. The in situ test was simulated by three-dimensional, discrete element software, and the relationship between fluid pressure and normal displacement is shown in Figure 1b.
The results show that the fracture width and permeability increases with an increasing water pressure in the fracture surface, while the fracture width and permeability decreases with decreasing water pressure. The change of water pressure will lead to the change of the fracture width in the process of fracture seepage. The change of the fracture width will also lead to a significant change in permeability, since the fracture surface permeability and aperture width have a cubic power relationship. Based on the verifying examples shown in the figure, the discrete element method can accurately simulate the physical mechanics phenomena in fracture seepage.

Hydraulic-Mechanical Coupling with Strength Reduction Method in DEM
Strength reduction method is typically applied in safety factor (F) calculations by progressively reducing the shear strength of the rock material to bring the slope to a state of limiting equilibrium (situation that actual strength state (τ a ) that corresponds to the shear strength (τ s )). An updated trial material parameter derived through the strength parameter was divided by the strength reduction coefficient F trial . The two strength reduction coefficients of the research object are in equilibrium state, and the non-equilibrium state can be obtained, respectively, by the reduction of the parameters. The safety factor was calculated using the "bracketing approach" (the least squares method), which is similar to the method proposed by Dawson (1999) [26]. The strength reduction coefficient at the time when the system in the limit equilibrium state and satisfy the deviation requirement is the safety factor. Thereby, the strength parameters C trial and Φ trial were calculated by: Besides, the way to determine the stability state of the research object during the strength reduction process by the unbalanced force ratio in Figure 2. The system in continuous plastic flow state means that there is the unbalanced force still in the system, causing the increasing displacement, which also makes it possible to get the safety factor by analyses of the changes of macroscopic displacement. The ratio of unbalanced force to the mean absolute value of force exerted by each surrounding zone, the unbalanced force is the net force acting on a grid point. Nr, characteristic response time, a representative number of steps (maximum limit of 50,000) that characterizes the response time of the system.

Dam Foundation Model
The Capulin San Pablo hydropower station is located on the Tarcoles River in Costa Rica. It mainly includes retaining structures, water discharge structures, a left bank water conveyance system, and a left bank ground workshop. The retaining structure is a roller-compacted-concrete gravity dam comprised of nine dam sections, with a crest elevation of 196.50 m, a maximum dam height of 48.5 m, a crest width of 7.4 m, and a dam axis length of 165.5 m. The exposed strata in the dam site is mainly composed of volcanic breccias as well quaternary residual slope gravel soil, collapse slope gravel, and stone with a small amount of clay. The lower slope of the riverbed is distributed with a gently inclined clayey tuff, the poor character of, which will affect the stability of the dam. In the foundation of the dam, weak structural planes, such as the slopes of the gentle dip faults (designated F11, F12, F2, and F3 and listed in Table 1) are widely developed. Specific rock mass parameters and dam structure parameters are shown in Tables 2 and 3, respectively. In the analysis of fractured rock mass flow problems or the stability of hydraulic structures, taking into account the seepage effect will be closer to the reality, moreover, it is of great significance considering the tortuosity and roughness of the fracture that improves seepage model. On this basis, we summarized the general process of analyzing these problems taking into account the hydraulic-mechanical coupling based on the strength reduction method as follows as Figure 2.

Dam Foundation Model
The Capulin San Pablo hydropower station is located on the Tarcoles River in Costa Rica. It mainly includes retaining structures, water discharge structures, a left bank water conveyance system, and a left bank ground workshop. The retaining structure is a roller-compacted-concrete gravity dam comprised of nine dam sections, with a crest elevation of 196.50 m, a maximum dam height of 48.5 m, a crest width of 7.4 m, and a dam axis length of 165.5 m. The exposed strata in the dam site is mainly composed of volcanic breccias as well quaternary residual slope gravel soil, collapse slope gravel, and stone with a small amount of clay. The lower slope of the riverbed is distributed with a gently inclined clayey tuff, the poor character of, which will affect the stability of the dam. In the foundation of the dam, weak structural planes, such as the slopes of the gentle dip faults (designated F11, F12, F2, and F3 and listed in Table 1) are widely developed. Specific rock mass parameters and dam structure parameters are shown in Tables 2 and 3, respectively.  A dam model was built using the discrete element method, as shown in Figure 3a, a Mohr-Coulomb constitutive model was adopted for the rock mass, and the structural plane was constructed using a Coulomb-slip constitutive model. The Capulin San Pablo hydropower station dam foundation is perpendicular to the river cross section (1-1 section) and the river cross section (2-2 section), as shown in Figure 3b. The main structural plane of the dam foundation and the segregation characteristics between the dam sections are shown in Figure 4. The model of the structural stability of the dam sections 1-3 is simulated in detail based on the developed structural plane of dam sections 1-3 by taking into account rock classes III1B, III 2B, IV2c, and V with other rock layers and structural faces C01, C05, C06, C10, and C18 with other structural effects. The monitoring points on the upper and lower dam and on both sides of the UG7 upper and lower disc layouts are shown in Figure 5.  Figure 4. The model of the structural stability of the dam sections 1-3 is simulated in detail based on the developed structural plane of dam sections 1-3 by taking into account rock classes III1B, III 2B, IV2c, and V with other rock layers and structural faces C01, C05, C06, C10, and C18 with other structural effects. The monitoring points on the upper and lower dam and on both sides of the UG7 upper and lower disc layouts are shown in Figure 5.

Hydraulic-Mechanical Coupling Calculation
After the preparation of initial material parameters, the mechanical, structural parameters of aimed research object, and the required fluid material parameters (including density, cohesion, bulk modulus, and joint hydraulic parameters, and the in situ stress and pore pressure), the model boundary pore pressure and pore pressure gradient were set as shown in Figure 6, namely P = ρgh.

Hydraulic-Mechanical Coupling Calculation
After the preparation of initial material parameters, the mechanical, structural parameters of aimed research object, and the required fluid material parameters (including density, cohesion, bulk modulus, and joint hydraulic parameters, and the in situ stress and pore pressure), the model boundary pore pressure and pore pressure gradient were set as shown in Figure 6, namely P = ρgh.
For the mechanical calculation, the model was calculated to equilibrium considering the weight and the water level force of the upstream and downstream dams (the surface force distribution is shown in Figure 6). For hydraulic calculations, the fracture seepage analysis of the dam foundation was carried out considering the influence of grouting. The fissure water pressure of the structural planes under the working condition of flood level was then obtained. The calculation proceeded until the model reached the equilibrium state, during which time the fracture surface water pressure on the structure remained at a fixed value. Then, the parameters of the dam and bedrock structural plane was reduced to obtain the safety factor under different working condition. The normal stress of the structural plane is considered in the model. The displacement of the monitoring points during the reduction process is recorded, and the safety factor of the dam is expressed according to the corresponding strength reduction coefficient when the displacement is suddenly changed.

Numerical Results Analysis about Dam Deformation
In flood level conditions, the monitoring point (dam body and points A and B of UG7) deformation curve with the strength reduction coefficient of the relationship is shown in Figure 7. Deformation characteristics shown in the figure indicate that the dam body and bedrock remain stable when the strength reduction factor is 2.80. However, the deformation of the monitoring point of the rock mass and the structural plane have dramatic changes when the strength reduction coefficient is 3.00. Therefore, it can be ascertained that the overall safety factor of the dam foundation For the mechanical calculation, the model was calculated to equilibrium considering the weight and the water level force of the upstream and downstream dams (the surface force distribution is shown in Figure 6). For hydraulic calculations, the fracture seepage analysis of the dam foundation was carried out considering the influence of grouting. The fissure water pressure of the structural planes under the working condition of flood level was then obtained. The calculation proceeded until the model reached the equilibrium state, during which time the fracture surface water pressure on the structure remained at a fixed value.
Then, the parameters of the dam and bedrock structural plane was reduced to obtain the safety factor under different working condition. The normal stress of the structural plane is considered in Energies 2017, 10, 1544 11 of 15 the model. The displacement of the monitoring points during the reduction process is recorded, and the safety factor of the dam is expressed according to the corresponding strength reduction coefficient when the displacement is suddenly changed.

Numerical Results Analysis about Dam Deformation
In flood level conditions, the monitoring point (dam body and points A and B of UG7) deformation curve with the strength reduction coefficient of the relationship is shown in Figure 7. Deformation characteristics shown in the figure indicate that the dam body and bedrock remain stable when the strength reduction factor is 2.80. However, the deformation of the monitoring point of the rock mass and the structural plane have dramatic changes when the strength reduction coefficient is 3.00. Therefore, it can be ascertained that the overall safety factor of the dam foundation is 3.00. It can also be seen from the calculation results that the deformation characteristics of monitoring points A and B in UG7 are basically consistent with the deformation characteristics of the dam monitoring points. The displacement distribution characteristics of the dam deformation under different strength reduction factors are shown in Figure 8a. It can be seen from the figure that, when the strength reduction coefficient is 2.00, 2.80, and 3.00, the deformations of the entire dam are in ranges 10-20, 50-60, and 90-105 mm, respectively. With an increasing reduction coefficient, the deformation characteristics of the rock mass near the top of dam sections 4-7 and UG7 first appear with the deformation's continuous expansion. The deformation of the entire dam foundation tends to be destroyed and destabilized from the lower level. It not only can be concluded that the shear failure is most likely to occur in dam sections 4-7 due to the influence of the UG7 clay tuff weak rock strata, but also that section 3 of the dam is easy to slide along the joint between sections 2 and 3.The dam foundation deformation and 1-1 dam base profile deformation distribution are shown in Figure 8b when strength reduction coefficient is 3.0.
The distribution of water pressure on the main structural surface under flood conditions is shown in Figure 9. Under the given seepage boundary condition, the distribution of fissure pressure in all structural planes is obtained after seepage calculation. The water head of the upstream side is higher and the fissure water pressure is 0.588 MPa.  The distribution of water pressure on the main structural surface under flood conditions is shown in Figure 9. Under the given seepage boundary condition, the distribution of fissure pressure in all structural planes is obtained after seepage calculation. The water head of the upstream side is higher and the fissure water pressure is 0.588 MPa.

Discussion
The Chinese standard entitled "Design Code for Concrete Gravity Dam [27]" stipulates that the safety factor needs to exceed 3.0 under basic load combination conditions when the strength of the rock mass and the structural plane is evaluated by the shear strength. In other countries, the design of anti-sliding stability of dam foundation is mainly referred to in Reclamation Bureau Report No. REC-ERC-74-10 [28]. This report states that the safety factor for normal load should be greater than 4.0, that the accident load should be greater than 2.7, and that the extreme load combination should be greater than 1.3.
Under flood conditions, the safety factor of the deep anti-sliding stability of the Capulin San Pablo dam basically meets the Chinese standard, but it does not meet the U. S. Bureau of Reclamation standard due to the fact that the security margin is very small. It is necessary to adopt more reliable measures, such as local replacement treatment, to increase the overall anti-sliding stability of the dam in order to meet the design standards in consideration of UG7, which is the most important factor affecting the stability of the dam foundation. In addition, the increase in uplift pressure of the dam foundation is caused by the increase of fissure water pressure in the structural

Discussion
The Chinese standard entitled "Design Code for Concrete Gravity Dam [27]" stipulates that the safety factor needs to exceed 3.0 under basic load combination conditions when the strength of the rock mass and the structural plane is evaluated by the shear strength. In other countries, the design of anti-sliding stability of dam foundation is mainly referred to in Reclamation Bureau Report No. REC-ERC-74-10 [28]. This report states that the safety factor for normal load should be greater than 4.0, that the accident load should be greater than 2.7, and that the extreme load combination should be greater than 1.3.
Under flood conditions, the safety factor of the deep anti-sliding stability of the Capulin San Pablo dam basically meets the Chinese standard, but it does not meet the U. S. Bureau of Reclamation standard due to the fact that the security margin is very small. It is necessary to adopt more reliable measures, such as local replacement treatment, to increase the overall anti-sliding stability of the dam in order to meet the design standards in consideration of UG7, which is the most important factor affecting the stability of the dam foundation. In addition, the increase in uplift pressure of the dam foundation is caused by the increase of fissure water pressure in the structural plane, which is the main factor affecting the stability of the dam foundation, and certain grouting measures may be considered to deal with it.
The traditional method of slope stability analysis is the limit analysis and the limit equilibrium method. The limit analysis method is based on the upper-and lower-bound theorems for the development of plasticity theory (Davis and Selvadurai, 2002) [29]. These theorems strictly limit the failure condition to a fully plastic material, and obeys the associated flow criteria, and a conservative estimate is thus given when calculating the safety factor. The limit equilibrium method needs to assume the shape and position of the slip surface. Although the limit equilibrium method adopts the upper-bound criterion of limit analysis, it is difficult to obtain the sliding surface of the rock and stable safety factor due to the fact that the upper-bound criterion's accuracy requirements are not satisfied. In the hybrid model, the discrete network model and the continuous medium model are used to simulate primary and secondary fracture, respectively. This model comes closer to simulating the real rock mass by combining the strengths of the two models. There is no doubt that the hybrid model is closer to reality under the condition of understanding the distribution of rock fissures. Wei andHudson (1993-1998) [30] established the hybrid model according to the discrete element and finite element methods in the near and far fields. Furthermore, they simulated the coupling behavior of the fractured rock mass, and the model is validated in the hydraulic-mechanical coupling analysis of the gravity dam foundation. However, further research is needed because of the difficulty of model building and the small sample size.
For the stability analysis of rock mass controlled by structural plane under hydraulic-mechanical interaction, the strength reduction method based on DEM in analysis of complex rock mass seepage is optimized as follows: (a) The discrete element method can fully reflect the deformation character caused by the difference between dam foundation rock parameters and structural plane parameters, which may affect the safety factor of local dam sections. (b) The strength reduction method in combination with the hydraulic-mechanical analysis produces an efficient solution that can express the continuous motion in the model, even in an unstable system. (c) The strength reduction method, based on the discrete element method, can automatically search for the minimum safety factor and the corresponding failure modes under the actual combination of structural planes, which is closer to a real-world situation.
A discrete fracture network can accurately reflect the seepage performance of a rock mass under the coupling of hydraulic-mechanical in fractured rock mass. However, when confronted with a large number of stochastic structures, the complicated process is difficult to evaluate. In this paper, it was assumed that the rock mass medium is impermeable and all of the fluid flows in the main control structure are suitable for the seepage model of the Capulin San Pablo dam. It is necessary to determine whether the hypothesis is satisfied if this method is used for stability analysis of other research objects.
The analysis of the stability of the dam based on the seepage of the fracture and using the strength reduction method can better and more comprehensively consider various factors. The safety factor of the dam obtained can be used to evaluate the anti-sliding stability of the dam and the stability analysis of the dam foundation slope. It can also provide the numerical simulation basis for effective and reasonable seepage control measures aimed at the weak parts of dam foundations.

Conclusions
A three-dimensional dam body model was built to reflect the dam shape and structures in dam and rock mass with discrete element software, designated the 3DEC model. The characteristics and applicable conditions of different seepage models were studied, and the safety factor of the dam was obtained using the strength reduction method. Finally, the stability of the dam foundation under seepage conditions was discussed. Several conclusions were reached, as follows: (a) The modified cubic law considered the tortuosity and roughness of the fracture that improves fracture network models, which provided an approach for accurately calculating hydraulic properties and made the analysis method of fractured rock mass more applicable. (b) The discrete element method was used to analyze the hydraulic-mechanical coupling problem in a fractured rock mass, and it was deemed suitable for analyzing such a problem. The results of analysis were more reasonable and were found to be closer to the actual problems in a fractured rock mass. (c) The strength reduction method based on the discrete element method has the advantage that it can be used to analyze the stability of the complex dam-foundation system. It can automatically search for the minimum safety factor and the failure mode under the condition of a combination of structural planes. This method can produce results more in line with the true state of dam foundation damage and is also suitable for the analysis of hydraulic-mechanical intersection failure under the action of a complex fractured system.