The Unsaturated Hydromechanical Coupling Model of Rock Slope Considering Rainfall Infiltration Using DDA

Water flow and hydromechanical coupling process in fractured rocks is more different from that in general porous media because of heterogeneous spatial fractures and possible fracture-dominated flow; a saturated-unsaturated hydromechanical couplingmodel using a discontinuous deformation analysis (DDA) similar to FEM and DEM was employed to analyze water movement in saturated-unsaturated deformed rocks, in which the Van-Genuchten model differently treated the rock and fractures permeable properties to describe the constitutive relationships. The calibrating results for the dam foundation indicated the validation and feasibility of the proposed model and are also in good agreement with the calculations based on DEM still demonstrating its superiority. And then, the rainfall infiltration in a reservoir rock slope was detailedly investigated to describe the water pressure on the fault surface and inside the rocks, displacement, and stress distribution under hydromechanical coupling conditions and uncoupling conditions. It was observed that greater rainfall intensity and longer rainfall time resulted in lower stability of the rock slope, and larger difference was very obvious between the hydromechanical coupling condition and uncoupling condition, demonstrating that rainfall intensity, rainfall time, and hydromechanical coupling effect had great influence on the saturatedunsaturated water flow behavior and mechanical response of the fractured rock slopes.


Introduction
Fractured rock masses including numerous discontinuities with various attitudes and different scales are complicated geological media that have undergone a long period of geological evolution.And the corresponding prediction and description of deformation, stress, and groundwater flow in the fractured rocks are of interest in many geotechnical areas such as slope stability [1,2] and seepage control for hydraulic projects [3,4].However, different from the general porous media, the permeability of intact rock is very low compared with that in fractures [5]; water flows predominantly through a single fracture or fracture networks [6]; it is more complex to describe the seepage behavior in such rocks.From the achievements about water flow in fractured rocks, numerical methods are often used to model seepage behaviors through fractured rocks, including the equivalent continuum approaches such as FDM, FEM [7], and BEM, the dual-continuum approach introduced by Barenblatt (1960) commonly used in mimicking water flow and mass transport in fractured porous medium [8,9], and discrete fracture network (DFN) [10][11][12][13].However, the continuum model is to view the fractured rock masses as an equivalent porous medium in which the heterogeneity of rock masses is averaged based on the concept of representative element volume (REV).Thus, it is desirable to evaluate spatial and random distribution of the fractured rock flow [12], and then a sound numerical method of intact rock and fractures should be proposed to model water movements in the deformed rocks for performance assessments and design optimization in a lot of rock engineering such as rock slopes and dam foundation.
And above researches are mainly about rock steady and nonsteady seepage focusing on water flow in the fully 2 Geofluids saturated mode; the unsaturated flow behavior in such fractured rock masses is very often encountered in geotechnical problems such as atomized rain or reservoir water level variation on the dam and foundation, embankment, and rainfall infiltrating into the rock slope.Whereas few investigations about saturated-unsaturated hydromechanical coupling model in the fractured rocks have been implemented, so the unsaturated flow behavior especially reflecting the hydromechanical coupling effect in fractured rocks cannot be ignored and should be deeply investigated considering the rock deformation.A discontinuous deformation analysis (DDA) [14] may be employed for time discretization and spatial approximation to solve the rock deformation.And in the DDA method, the displacements and strains of the blocks are determined by solving a system of simultaneous equilibrium equations which are formulated by minimizing the total potential energy of the block system, similar to the conventional FEM.In addition, the DDA method allows the blocks to have various modes of motion, such as sliding, rotation, and opening, so that large displacements can be simulated, which is similar to the DEM [15,16].Meanwhile, with continuous contributions to its improvements and applications, the DDA was well developed in terms of both theory and computing codes [17][18][19][20][21] to analyze the mechanical response of the rocks with random fractures or faults.And then, in order to determine the effective saturation and relative permeability [22,23], the Van-Genuchten model was written in saturated-unsaturated seepage codes differently treating the rock and fractures permeable properties describing the constitutive relationships during the unsaturated flow processes.Thereby, combined with the DDA and seepage codes of fractured rocks, a saturated-unsaturated hydromechanical coupled model using DDA has been proposed.
For the sake of better applying the proposed coupled model to the corresponding practical engineering, the validity of the developed method was verified by the consistency of numerical results and computing results by the other method DEM in the analysis of a dam foundation.And also, a fractured rock slope considering rainfall infiltration as a case study and the water pressure distribution on the fault surface and inside the rocks was demonstrated with rainfall time increasing, and also the displacement, principle stress distribution, and slope stability were described considering the hydromechanical coupling and uncoupling process, indicating that rainfall intensity, rainfall time, and hydromechanical coupling effect under unsaturated conditions have great influence on the movement of seepage flow through fractured rock slopes and corresponding stability.

Unsaturated Hydromechanical Coupling Model of Fractured Rock Mass Using DDA
In view of the complexity of fractured network and hydraulic parameters determination, some traditional analysis models are limited for saturated-unsaturated fluid flow of fractured rocks.For better describing the hydraulic and mechanical characteristic of the rock masses, an unsaturated hydromechanical coupling model based on DDA was brought forward to simulate the water pressure, fracture deformation, and stress distribution to analyze the seepage field and stress field variation.

Unsaturated Flow Model of Fractured Rock Masses.
It is known that the entire area of the fractured rocks can be divided into saturated domains and unsaturated domains.The movement of partially saturated flow in an unsaturated domain appears because of the capillary effect.The governing equation of unsaturated water flow in porous media can be written as where  is water density; S is water resource;  is porosity;   is saturation, 0 ≤   ≤ 1.0;  is water content,   = ; V  is Darcy velocity.
According to Darcy's law, where () is permeability tensor related to water content ; ℎ is head;   is the saturated permeability tensor of fractured rocks (,  = 1, 2, 3).Substituting ( 2) into (1) gives the unsaturated flow expressions in the fractured rocks as follows: where   is the relative permeability, 0 ⩽   < 1 in unsaturated field, and   = 1 in saturated domains; ℎ  is pressure head and ℎ =  + ℎ  ; (ℎ  ) is water capacity and (ℎ  ) = 0 in saturated domains;  is special parameter,  = 0 in unsaturated domains, and  = 1 in saturated filed;   is the storage coefficient,   = 0 for the unsaturated rock, and   is constant for saturated rocks.
To solve (3) in the saturated-unsaturated flow analysis, initial and boundary conditions are needed as follows.
(1) Initial condition is expressed as (

2.2.
The Theory of DDA Method.Discontinuous deformation analysis (DDA) is a numerical method that can be used to simulate the discontinuous deformation behavior of jointed rocks [14].The method not only can model all kinds of blocks cut by the fractures, but also deals with the contact [24] and restraint between all the blocks.And the method determined by minimizing the total potential energy [19] of all blocks is similar to FEM and allowing all the blocks to have various mechanical motions is similar to DEM.

Motion Equation of Block
System.Suppose  blocks in the block system, the potential energy of a block, includes the elastic strain energy Π  , initial stress potential energy Π  , point force potential energy Π  , linear force potential energy Π  , body force potential energy Π V , inertial energy Π  , and displacement constraint energy Π  can be written as According to the principle of minimum potential energy, derive (7) to obtain the stiffness matrix expression [  ] and load matrix [  ] of a single block, where  = 1, . . ., 6 and  = 1, . . ., 6. Thus, the equations of motion of the blocks can be given in the matrix form,

Treatment of Water Load.
Water flow behavior in the fractured rocks is significant for building the unsaturated hydromechanical coupling model; water pressure should be correctly determined on the blocks.Shown in Figure 1 representing the water pressure on a block, suppose  1 = { 1 ,  1 } and  2 = { 2 ,  2 } at points ( 1 ,  1 ), ( 2 ,  2 ), respectively, and the water pressure   = {  (),   ()} at any point (, ) can be described by (10), where 0 ≤  ≤ 1 Thus, the energy Π  of th block considering the water pressure as the linear force can be expressed by Then, derive (11) to get the right load item   , Assuming the constant constraint element, the load item   can be written as To improve the computing precision, the higher-order nonconstant constraint element is introduced to determine the load items with the following two ways.
(1) Assuming the regular rectangular element of a block with four points (  ,   ),  = 1, 2, 3, 4 under global coordinate; two transformation matrices are written as ( 14) and (15).And Substituting ( 15) into (12) gives the corresponding load item   as ( 16) (2) In general, not all the block elements are regular so that complex nonregular element is needed to be deeply developed.The displacement function will be expressed by isoparametric element function; the load item   can be written by where   (, ),  = 1, 2, 3, 4 is weight function.Suppose global coordinate of (, ) and local coordinate of (, ) and the relationship of (, ) in (, ) can be described using the shape function   (, ),  = 1, 2, 3, 4, written as where   ,   ,  = 1, 2, 3, 4, are the vertex coordinates of a block element and {, }  can be expressed by (20).And substituting ( 21) into (19) gives the expressions of weight function   (, ),  = 1, 2, 3, 4, written as ( 22) and shape function [fun] 12×8 as (23); the load item   can be expressed by (24).The expression of the water load on the block system using DDA can realize the hydromechanical coupling effect coming from the water flow inside the block and the block deformation where

Unsaturated Hydromechanical Coupling Model of Fractured Rock Masses
Using DDA.The unsaturated coupled model using DDA has been proposed to solve the water pressure variation and mechanical response in saturatedunsaturated rocks.Note that the aperture variation of rock fractures resulting in obvious water pressure variation and rock deformation, corresponding real-time coupling process cannot be ignored.However, the real fractures are coarse and nonparallel, Cubic Law cannot be directly applied in this study.For the sake of getting new permeability tensor at any time, fractures deformation including normal deformation and shear deformation [24] should be completely dealt based on DDA calculation.As mentioned in previous studies about water flow in the fractured rocks, natural fractures cannot completely close under high normal stress with enough mechanical aperture.Assuming  01 ,  02 are the initial fractures aperture for each side of the fracture,  1 ,  2 representing the fracture aperture after rock deformed are written by (25) considering fracture normal deformations  1 and  2 (open is positive) and shear deformations  1 and  2 computed by DDA.And then, the corresponding flow rate can be expressed by (26); thus the equivalent hydraulic aperture can be calculated by ( 27) to get new permeability tensor at any time where   is the dilation angle;  =  1 / 2 ;  is the flow rate;   is the flow rate corresponding to average aperture   ; and Therefore, the proposed coupled model can consider the block normal deformation and shear dilation when water is flowing through the rock fractures and can also be applied to determine the unsaturated hydromechanical coupling behavior in the fractured rocks.Shown in Figure 2 is the flow chart of the proposed coupled model using DDA, recording the whole calculating process.

Calibration of the Proposed Coupled Model Based on DDA
For the sake of calibrating the proposed hydromechanical coupling model, shown in Figure 3  dam heel as the original of the coordinates with -axis directing toward the downstream and -axis pointing upwards.The mechanical parameters of dam are the density 2400 kg/m 3 , Young's module 19.6 GPa, and Poisson's ratio 0.167.And corresponding mechanical parameters of rock foundation are the density 2650 kg/m 3 , Young's module 8.0 GPa, Poisson's ratio 0.24, internal friction angle 45 ∘ , and cohesion 2 MPa.Suppose a group of orthogonal fractures in the rock foundation with aperture of 1 mm, the normal stiffness, and tangential stiffness of 20 GPa/m.The boundary conditions include the horizontal displacement restraint for the upstream and downstream boundaries of the rock foundation and vertical displacement restraint for the bottom boundary of the rock foundation, and water level acting on the upstream face of the dam body is 50 m and on the downstream face is 0 m.And then, both of the numerical models based on DDA and DEM methods were employed to make hydromechanical coupling calculations considering the horizontal water pressure, uplift pressure, and water flow force in rock foundation.For the sake of comparison, Figures 3, 4, and 5 give the calculating results of DDA and DEM models, respectively.The water head (unit: m), displacement (unit: m) and maximum principal stress (unit: Pa) based on proposed DDA model are shown in Figures 4(a), 5(a), and 6(a), indicating that the fractures in the upstream are tending to be tensile causing the fractures permeability increasing; however, the fractures in the downstream are tending to be compressed causing the fractures permeability decreasing.Therefore, the curves of water head in upstream move toward the upstream, and the curves of water head in downstream move toward the downstream; the variations are consistent with real seepage characteristics in the rocks.And it can be seen from the computing results shown in Figures 4(b), 5(b), and 6(b) based on DEM that the variations are nearly identical as determined by DDA.The above comparison of two methods demonstrated the feasibility and reliability of the proposed coupled model using DDA and can also be applied to analyze the mechanical characteristics of rock masses.In addition, the proposed model considering the heterogeneous deformation of the fracture deformation with different permeability and the better contact criterion and convergence criterion [19] will get more accurate computing results close to the real performance, so the model is superior to DEM to solve the practical hydromechanical problems in geotechnical engineering.

Analysis of a Fractured Rock Slope Considering Rainfall Infiltration
As noted in previous investigations of unsaturated fluid flow in slope engineering [1,2,25], rainfall infiltration will increase the water pressure inside the rocks and fractures and weaken the rock strength parameters; the landslide often occurs in the rain season.And the analysis of the fractured rock slope has been employed to demonstrate the capacity and robustness of the proposed coupled model for describing the rainfall infiltrating process influencing the slope water pressure, displacement, and stress distribution.Different from the saturated steady seepage or unsaturated in porous media, rainfall infiltrating into the rocks will flow through the fractures, directly loading on the fractures surface causing the real-time variation of the hydraulic aperture; the permeability of the rock and fractures will change accordingly, which also implies the unsaturated hydromechanical coupling process.As shown in Figure 7, a rock slope with a height of 180.0 m, length of 20.0 m, and width of 1.0 m is considered.There is a set of obvious faults with initial hydraulic aperture of 1 mm developed in the rock slope.A series of computing results about water pressure inside the rock and on the fault surface, displacement, stress, and stability coefficient under different loading conditions have been analyzed in detail based on the proposed DDA method.

Case Introduction.
As an initial condition, the rock slope is treated as a water head boundary assuming an external water level 110.0 m on the upstream face and 32.0 m on the downstream face, and the bottom boundaries were impermeable.The rainfall infiltration boundaries are imposed on the slope surface.And the initial saturated permeability coefficient of the rocks is assumed to be 1.39 × 10 −5 m/s.For the unsaturated coupled model, water retention curves are very important to express the unsaturated rock characteristics; particularly the Van-Genuchten model (V-G) is often used to describe the constitutive relations such as -  and -ℎ.Table 1 gives the constitutive relations based on V-G model [26].And Table 2 gives the mechanical parameters of the rock and the fault.figures that, due to the rainfall infiltration, water content of the rock slope is increasing and the matric suction is reducing gradually, and the negative pressure region continues to reduce progressively; the expansion of the transient area leads to large increase of water pressure inside the rock slope, which can drive the wetting front quickly into a deeper location.

Analysis for the
In particular, with rainfall infiltrating inside the slope, there will appear greater positive transient water pressure and larger range of transient saturated zone on the rock slope surface.The results also show that transient saturated zone will quickly be larger in shorter period of time with greater rainfall intensity.In addition, Figure 10 presents the water pressure (unit: m) distribution on the fault surface at different rainfall time, indicating the flow behaviors of rainfall infiltration inside the fault.It can be observed from Figure 10(a) with rainfall intensity of 8 mm/h that water pressure on the fault surface varies slightly, whereas it varies more greatly with rainfall intensity of 20 mm/h shown in Figure 10(b) considering the same rainfall infiltrating time.Meanwhile, as for the larger range of the unsaturated zone inside the rock slope, greater rainfall intensity will cause more and more water infiltration into the unsaturated zone to become the saturated zone, and the groundwater level will ascend quickly so as to increase the water pressure on the fault surface.Therefore, the water pressure on the sliding surface in the unsaturated zone will increase when rain falling time lasts with same rainfall intensity.And it can be seen from the two figures that, in the same period of rainfall infiltrating time, greater rainfall intensity will cause larger water content of the rocks in a shorter period resulting in greater unit weight increasing, and the water pressure head will rise more obviously so that the water pressure in the transient unsaturated zone increases greater and more quickly, which demonstrates that the rock slope will be sliding more easily considering rainfall infiltration of greater rainfall intensity.
From the above calculating results in agreement with the water flow behaviors in the saturated-unsaturated rocks, it can be noticed that the proposed coupled model is also applicable to analyze the unsaturated hydromechanical coupling characteristics under complex conditions.

Comparison under Unsaturated Hydromechanical Coupling and Uncoupling Conditions.
Based on the deeper calibration of the proposed coupled model in Section 4.2, further researches will be advanced to clearly compare the mechanical characteristic of the rock slope under unsaturated hydromechanical coupling conditions and uncoupling conditions.Figures 11 and 12 show the water pressure head distribution under coupling conditions and uncoupling conditions at time T = 60 h and T = 80 h with the rainfall intensity of 8 mm/h and 20 mm/h.It can be seen from the above figures that the water pressure on the fault surface will be larger in the transient saturated zone considering the coupling process.And the comparison under two conditions demonstrates that coupling process causes water pressure head increase more obviously with longer rainfall infiltrating time and greater rainfall intensity, which is more unfavorable for the slope stability.
Furthermore, shown in Figures 13 and 14 are the maximum principle stress contours (unit: MPa) and displacement contours (unit: m) under coupling conditions and uncoupling condition considering the rainfall intensity of 20 mm/h at time  = 20 h,  = 60 h, and  = 80 h.It can be noted from the comparison of the figures that the  stress and displacement of the rock slope are larger under coupling condition than that under uncoupling conditions, which also indicates that the unsaturated hydromechanical effect will enlarge failure probability of the slope and should not be ignored especially under rainfall infiltration conditions.
For the sake of explaining the coupling effect influencing on the rock slope stability, Table 3 gives the stability coefficient of the rock slope under different conditions.The results indicate that the safety coefficients under hydromechanical coupling conditions are smaller than that under uncoupling conditions, which is in agreement with the water pressure, stress, and displacement variation.It can be concluded that the rainfall intensity, rainfall time, and hydromechanical coupling effect have significant influence on the slope stability.
Therefore, the hydromechanical coupling process under rainfall infiltration for a rock slope will cause the transient saturated zone to be larger, and the transient water pressure will accordingly increase the downslide force and water pressure causing the effective normal stress on the sliding surface reducing too much, which is more unfavorable for the rock slope stability.The comparison between uncoupling conditions and coupling conditions proves the applicability of the unsaturated hydromechanical coupled model based DDA to analyze the mechanical response of related practical projects.

Conclusions
DDA similar to FEM and DEM has attracted increasingly wide attention in the geotechnical engineering; it can be applied to analyze the rock deformation and fractures aperture variation.And for the hydromechanical coupling issues of the fractured rocks, the fractures aperture variation changing the permeability will seriously influence the water pressure distribution; DDA method is more favorable to solve the fractures deformation under complex conditions.Therefore, an unsaturated coupled model based on DDA has been proposed for the analysis of saturated-unsaturated water flow behavior and mechanical response of the fractured rocks such as dam foundation and rainfall infiltration in rock slopes, and a systematic numerical code was developed to analyze the hydromechanical process, indicating the feasibility and

Figure 1 :
Figure 1: Water pressure on the block.
and   are normal flow rate;   is direction cosine of external normal;  0 is initial time; Γ 1 is head boundary; Γ 2 is flow rate boundary; Γ 3 is saturated overflow boundary; Γ 4 is unsaturated overflow boundary.

Table 1 :
Seepage parameters of rocks.

Table 2 :
Mechanical parameters of rocks.

Table 3 :
Comparison of factors of stability for different conditions.