Table of Contents Author Guidelines Submit a Manuscript
Volume 2017, Article ID 1513421, 15 pages
Research Article

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

1Key Laboratory of New Technology for Construction of Cities in Mountain Area, Chongqing University, Ministry of Education, Chongqing 400045, China
2School of Civil Engineering, Chongqing University, Chongqing 400045, China

Correspondence should be addressed to Xianshan Liu; moc.361@uomuomzl

Received 21 September 2017; Revised 14 November 2017; Accepted 22 November 2017; Published 19 December 2017

Academic Editor: Qinghui Jiang

Copyright © 2017 Xianshan Liu and Ming Xu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 coupling model 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 saturated-unsaturated water flow behavior and mechanical response of the fractured rock slopes.

1. 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) [1013]. 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 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 [1721] 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.

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

2.1. 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, ; is water content, ; 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 .

Substituting (2) into (1) gives the unsaturated flow expressions in the fractured rocks as follows: where is the relative permeability, in unsaturated field, and in saturated domains; is pressure head and ; is water capacity and in saturated domains; is special parameter, in unsaturated domains, and in saturated filed; is the storage coefficient, 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.

Initial condition is expressed as

The boundary conditions generally including the water head, flux, and seepage-face boundaries are described by where and are normal flow rate; is direction cosine of external normal; is initial time; is head boundary; is flow rate boundary; is saturated overflow boundary; is unsaturated overflow boundary.

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.

2.2.1. Displacement Functions of Deformable Blocks

Different from FEM selecting the node displacements as the degrees of freedom, DDA selects the rigid body displacements and the strains of block elements as the basic unknown variables to solve for solutions. DDA conducts calculations based upon time steps in which the large displacement and deformation are the accumulating results of the small ones obtained from time-stepwise calculations. Suppose that an arbitrarily shaped polygonal block has uniform stress and strain, the unknown variables including parallel movement , rotation , normal strain, and tangential strain of the block can be written as where are, respectively, shape matrix and displacement matrix.

2.2.2. 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 , 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 and

Thus, the equations of motion of the blocks can be given in the matrix form,

2.2.3. 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 and at points , respectively, and the water pressure at any point can be described by (10), where

Figure 1: Water pressure on the block.

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.

Assuming the regular rectangular element of a block with four points under global coordinate; two transformation matrices are written as (14) and (15). And Substituting (15) into (12) gives the corresponding load item as (16)

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 is weight function. Suppose global coordinate of and local coordinate of and the relationship of in can be described using the shape function , written as where , , 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 , written as (22) and shape function 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 and

2.3. 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 saturated-unsaturated 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 , are the initial fractures aperture for each side of the fracture, , representing the fracture aperture after rock deformed are written by (25) considering fracture normal deformations and (open is positive) and shear deformations and 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; ; 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.

Figure 2: Flow chart of the unsaturated hydromechanical coupling analysis using DDA.

3. Calibration of the Proposed Coupled Model Based on DDA

For the sake of calibrating the proposed hydromechanical coupling model, shown in Figure 3 is a model of concrete gravity dam with rock foundation that measures 50 m in height, 6 m in width at dam top, and 30 m in width at dam bottom, with the upstream face vertical and the downstream one having a slope of 1 : 0.682. The rock foundation is 105 m in length and 50 m in height. The calculation domain takes the 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/m3, Young’s module 19.6 GPa, and Poisson’s ratio 0.167. And corresponding mechanical parameters of rock foundation are the density 2650 kg/m3, 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.

Figure 3: Block element meshes for DDA analysis (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.

Figure 4: Water head contours for coupling (solid line) and initial condition (dashed line) (m).
Figure 5: Vertical displacement contours under coupling conditions.
Figure 6: Maximum principle stress contours under coupling conditions.

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

Figure 7: Diagram of a rock slope (unit: m).
4.1. 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.

Table 1: Seepage parameters of rocks.
Table 2: Mechanical parameters of rocks.
4.2. Analysis for the Rock Slope Considering Rainfall Infiltration

Suppose that the rainfall lasts for 80 h with rainfall intensity of 8 mm/h and 20 mm/h, and analyze the water pressure (unit: m) at different time. The distribution of water pressure at different time is shown in Figures 8 and 9, respectively. Negative water pressure indicates unsaturated area; positive water pressure demonstrates fully saturated area, which is in accordance with the rainfall filtration mechanism of the rock slope. It can be observed from these 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. 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.

Figure 8: Pressure head contours of rainfall intensity 8 mm/h.
Figure 9: Pressure head contours of rainfall intensity 20 mm/h.

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.

Figure 10: Pressure head contours on fault considering rain falling time.

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.

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

Figure 11: Pressure head for coupling (dashed line) and uncoupling (solid line) (8 mm/h).
Figure 12: Pressure head contours for coupling (dashed line) and uncoupling (solid line) (20 mm/h).

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.

Figure 13: Maximum principle stress contours for coupling (dashed line) and uncoupling (solid line) (20 mm/h).
Figure 14: Total displacement contours for coupling (dashed line) and uncoupling (solid line) (20 mm/h).

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.

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

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.

5. 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 applicability of the proposed coupled model. On the basis of this work, the following conclusions can be drawn:

For better modelling the deformation of the rocks and fractures, the higher-order isoparametric element is selected in DDA procedure to improve the numerical precision. And on the basis of Cubic Law, the fracture aperture under complex conditions can be derived considering the fracture norm deformation and shear deformation to calculate the permeability at any time. And also selecting the V-G model for describing the constitutive relations and appropriate evaluation of the capacity term, the unsaturated coupled model has been proposed to describe the hydromechanical response of the fractured rocks. Meanwhile, the effectiveness of the model applying in a dam foundation has been calibrated by the satisfactory consistency of the results by DDA and computing results by DEM.

The proposed model has been applied to describe the water variation at various time under a rainfall infiltration process in a rock slope and the unsaturated flow behavior with respect to the variation of water pressure inside the rocks and on the fault surface is well revealed. It can also be noted that the matric suction of the rock slope is smaller and smaller when rainfall infiltrating lasts, and greater rainfall intensity causes safety coefficient to reduce quickly, which is in good agreement indicating the feasibility of the proposed model.

Further research has been advanced by comparison of the pressure head, principle stress, and displacement under coupling conditions and uncoupling conditions, demonstrating that the above corresponding values are increasing and the safety coefficients are reducing at the same time under coupling conditions, and also indicating that the rainfall intensity, rainfall time, and hydromechanical coupling effect under unsaturated conditions have great influence on the stability of the rock slope. It also explains that it is more dangerous under rainfall infiltration in a rock slope considering the hydromechanical coupling process. Therefore, the proposed unsaturated coupled model based on DDA can better describe the hydromechanical response and the stability of the rock slope.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The research is financially supported by the National Natural Science Foundation of China nos. 51779021 and 51478065 and Fund of National Engineering and Research Center for Mountainous Highways no. GSGZJ-2015-09.


  1. F. Cai and K. Ugai, “Numerical analysis of rainfall effects on slope stability,” International Journal of Geomechanics, vol. 4, no. 2, pp. 69–78, 2004. View at Publisher · View at Google Scholar · View at Scopus
  2. S. E. Cho and S. R. Lee, “Instability of unsaturated soil slopes due to infiltration,” Computers & Geotechnics, vol. 28, pp. 185–208, 2001. View at Publisher · View at Google Scholar
  3. Y. Chen, R. Hu, W. Lu, D. Li, and C. Zhou, “Modeling coupled processes of non-steady seepage flow and non-linear deformation for a concrete-faced rockfill dam,” Computers & Structures, vol. 89, no. 13-14, pp. 1333–1351, 2011. View at Publisher · View at Google Scholar · View at Scopus
  4. Y. Chen, R. Hu, C. Zhou, D. Li, and G. Rong, “A new parabolic variational inequality formulation of Signorini's condition for non-steady seepage problems with complex seepage control systems,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 35, no. 9, pp. 1034–1058, 2011. View at Publisher · View at Google Scholar · View at Scopus
  5. K. Pruess, “On water seepage and fast preferential flow in heterogeneous, unsaturated rock fractures,” Journal of Contaminant Hydrology, vol. 30, no. 3-4, pp. 333–362, 1998. View at Publisher · View at Google Scholar · View at Scopus
  6. L. J. Pyrak-Nolte, C. D. Montemagno, and D. D. Nolte, “Volumetric imaging of aperture distributions in connected fracture networks,” Geophysical Research Letters, vol. 24, no. 18, pp. 2343–2346, 1997. View at Publisher · View at Google Scholar · View at Scopus
  7. G. Gottardi and M. Venutelli, “UPF: Two-dimensional finite-element groundwater flow model for saturated-unsaturated soils,” Computers & Geosciences, vol. 27, no. 2, pp. 179–189, 2001. View at Publisher · View at Google Scholar · View at Scopus
  8. T. Graf and R. Therrien, “Variable-density groundwater flow and solute transport in irregular 2D fracture networks,” Advances in Water Resources, vol. 30, no. 3, pp. 455–468, 2007. View at Publisher · View at Google Scholar · View at Scopus
  9. C. Masciopinto and M. C. Caputo, “Modeling unsaturated-saturated flow and nickel transport in fractured rocks,” Vadose Zone Journal, vol. 10, no. 3, pp. 1045–1057, 2011. View at Publisher · View at Google Scholar · View at Scopus
  10. Q. Jiang, C. Yao, Z. Ye, and C. Zhou, “Seepage flow with free surface in fracture networks,” Water Resources Research, vol. 49, no. 1, pp. 176–186, 2013. View at Publisher · View at Google Scholar · View at Scopus
  11. C. Yao, Q. Jiang, W. Wei, and C. Zhou, “The variational inequality formulation for unconfined seepage through three-dimensional dense fracture networks,” Science China Technological Sciences, vol. 56, no. 5, pp. 1241–1247, 2013. View at Publisher · View at Google Scholar · View at Scopus
  12. Z. Ye, Q. Jiang, C. Zhou, and Y. Liu, “Numerical analysis of unsaturated seepage flow in two-dimensional fracture networks,” International Journal of Geomechanics, vol. 17, no. 5, Article ID 04016118, 2017. View at Publisher · View at Google Scholar · View at Scopus
  13. L. Jing and J. A. Hudson, “Numerical methods in rock mechanics,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 4, pp. 409–427, 2002. View at Publisher · View at Google Scholar · View at Scopus
  14. G. H. Shi, Discontinuous deformation analysis-a new numerical model for the statics and dynamics of block system, University of California, Berkeley, Calif, USA, 1988.
  15. J.-H. Wu, W.-K. Lin, and H.-T. Hu, “Assessing the impacts of a large slope failure using 3DEC: The Chiu-fen-erh-shan residual slope,” Computers & Geosciences, vol. 88, pp. 32–45, 2017. View at Publisher · View at Google Scholar · View at Scopus
  16. E. G. Nezami, Y. M. A. Hashash, D. Zhao, and J. Ghaboussi, “A fast contact detection algorithm for 3-D discrete element method,” Computers & Geosciences, vol. 31, no. 7, pp. 575–587, 2004. View at Publisher · View at Google Scholar · View at Scopus
  17. Y. Ning, J. Yang, G. Ma, and P. Chen, “Modelling rock blasting considering explosion gas penetration using discontinuous deformation analysis,” Rock Mechanics & Rock Engineering, vol. 44, no. 4, pp. 483–490, 2011. View at Publisher · View at Google Scholar · View at Scopus
  18. Y. Jiao, X. Zhang, and J. Zhao, “Two-dimensional DDA contact constitutive model for simulating rock fragmentation,” Journal of Engineering Mechanics, vol. 138, no. 2, pp. 199–209, 2012. View at Publisher · View at Google Scholar · View at Scopus
  19. Q. Jiang, Y. Chen, C. Zhou, and M.-C. R. Yeung, “Kinetic energy dissipation and convergence criterion of discontinuous deformations analysis (DDA) for geotechnical engineering,” Rock Mechanics and Rock Engineering, vol. 46, no. 6, pp. 1443–1460, 2013. View at Publisher · View at Google Scholar · View at Scopus
  20. Q. Jiang and C. Zhou, “A rigorous solution for the stability of polyhedral rock blocks,” Computers & Geosciences, vol. 90, pp. 190–201, 2017. View at Publisher · View at Google Scholar
  21. N. Chen, J. Kemeny, Q. Jiang, and Z. Pan, “Automatic extraction of blocks from 3D point clouds of fractured rock,” Computers & Geosciences, vol. 109, pp. 149–161, 2017. View at Publisher · View at Google Scholar
  22. R. Hu, Y.-F. Chen, H.-H. Liu, and C.-B. Zhou, “A water retention curve and unsaturated hydraulic conductivity model for deformable soils: Consideration of the change in pore-size distribution,” Géotechnique, vol. 63, no. 16, pp. 1389–1405, 2013. View at Publisher · View at Google Scholar · View at Scopus
  23. R. Hu, Y. Chen, H. Liu, and C. Zhou, “A relative permeability model for deformable soils and its impact on coupled unsaturated flow and elasto-plastic deformation processes,” Science China Technological Sciences, vol. 58, no. 11, pp. 1971–1982, 2015. View at Publisher · View at Google Scholar · View at Scopus
  24. A. R. Keneti, A. Jafari, and J.-H. Wu, “A new algorithm to identify contact patterns between convex blocks for three-dimensional discontinuous deformation analysis,” Computers & Geosciences, vol. 35, no. 5, pp. 746–759, 2008. View at Publisher · View at Google Scholar · View at Scopus
  25. P. H. S. W. Kulatilake, L. Wang, H. Tang, and Y. Liang, “Evaluation of rock slope stability for yujian river dam site by kinematic and block theory analyses,” Computers & Geosciences, vol. 38, no. 6, pp. 846–860, 2011. View at Publisher · View at Google Scholar · View at Scopus
  26. K.-H. Li and J.-R. Chai, “Research on capillary pressure-saturation curve for a single fracture in an unsaturated fractured mass,” Rock and Soil Mechanics, vol. 7, no. 8, pp. 1253–1257, 2006. View at Google Scholar · View at Scopus