# The Parabolic Variational Inequalities for Variably Saturated Water Flow in Heterogeneous Fracture Networks

Authors: Zuyang Ye; Qinghui Jiang; Chi Yao; Yanzhang Liu; Aiping Cheng; Shibing Huang; Yaning Liu
Journal: Geofluids
 (2018)

---

## Abstract

Fractures are ubiquitous in geological formations and have a substantial influence on water seepage flow in unsaturated fractured rocks. While the matrix permeability is small enough to be ignored during the partially saturated flow process, water seepage in heterogeneous fracture systems may occur in a non-volume-average manner as distinguished from a macroscale continuum model. This paper presents a systematic numerical method which aims to provide a better understanding of the effect of fracture distribution on the water seepage behavior in such media. Based on the partial differential equation (PDE) formulations with a Signorini-type complementary condition on the variably saturated water flow in heterogeneous fracture networks, the equivalent parabolic variational inequality (PVI) formulations are proposed and the related numerical algorithm in the context of the finite element scheme is established. With the application to the continuum porous media, the results of the numerical simulation for one-dimensional infiltration fracture are compared to the analytical solutions and good agreements are obtained. From the application to intricate fracture systems, it is found that water seepage flow can move rapidly along preferential pathways in a nonuniform fashion and the variably saturated seepage behavior is intimately related to the geometrical characteristics orientation of fractures.


---

## Body

## 1. Introduction

In the past decades, the understanding of water seepage flow in unsaturated, fractured rocks has been investigated by many researchers. Wang and Narasimhan [1] developed a conceptual model to simulate partially saturated flow in fractured porous media; the constitutive relationships and contact area were both summarized as functions of pressure head on the basis of a gamma aperture distribution. Peters and Klavetter [2] proposed a continuum model to evaluate the water movement in an unsaturated fractured porous medium, in which unsaturated hydraulic properties were obtained by a macroscopic approach and a microscopic approach. In order to describe the infiltration of a liquid front in an unsaturated fractured porous medium, an analytical solution was obtained by Nitao and Buscheck [3]. Therrien and Sudicky [4] developed a discrete fracture, saturated-unsaturated numerical model to analyze variably saturated flow and solute transport in fractured porous media, where the matrix and fractures were treated as three-dimensional and two-dimensional geometries, respectively. Other studies on variably saturated flow include those of Abdel-Salam and Chrysikopoulos [5], who used a mathematical model to investigate unsaturated flow in a quasi-three-dimensional fractured rock matrix system, and Masciopinto and Caputo [6], who suggested a numerical model for variably saturated flow and nickel transport through fractures by coupling plastic ring infiltrometer tests with time-lapse electrical resistivity tomography.In these conventional models, the macroscale continuum and dual-continuum concepts [2, 7] are frequently employed to describe water seepage flow in saturated-unsaturated fractured media. Under unsaturated conditions, the hydraulic behavior for heterogeneous fracture and matrix is homogenized in the manner of smooth sheets with very low flow rate and strong capillary suction effects into the partially saturated matrix that is performed based on the volume-averaging concepts [1, 3].Furthermore, it has been revealed by mounting evidence that the presence of fractures in unsaturated fractured rock formations can enhance the permeability of rock masses and that preferential and fast flow is associated with fractures [8–11]. From the previous studies [1, 12], matrix permeability of intact rocks such as welded tuffs is extremely low and generally several orders of magnitude smaller than that of fractures. According to a number of numerical calculations that consider water imbibition into wall rocks, the seepage flow over short time periods almost cannot be affected by rock matrix permeability in most simulations. Thus, the geometrical feature of the fracture networks can have a substantial influence on the variably saturated seepage flow and fluid transport, which highly depend on the connectivity and spatial distribution of fracture aperture, extent, and orientation. Additionally, the seepage and infiltration boundaries in the fracture network domain may be highly nonlinear, which makes it more difficult to develop theoretical and numerical solutions associated with the unsaturated water flow movement (Hu et al., 2016)[13, 14]. Recently, the PVI method has been extensively used on unsteady seepage problems in fracture networks and unsaturated flow problems in porous media [15] (Hu et al., 2016), where the boundary conditions on the seepage face and soil-atmosphere interfaces are described by a generalized complementary function. In order to model the water movement in unsaturated fracture networks, Ye et al. [16] proposed a finite element procedure for Richards’ equation and the related partial differential equations on the boundary conditions. Even though the Signorini-type equations on the seepage face or the infiltration surface were formulated by Ye et al. [16], the seepage points or the critical value between rainfall intensity and ponding depth is singular in mathematics, which greatly increases the nonlinearity in modeling the variably saturated flow problem. This issue has not been solved theoretically yet.In this study, the PVI method is extended by developing a systematic approach on modeling variably saturated water flow in the heterogeneous fracture networks and eliminating the singularity on the Signorini-type conditions, where the equivalent PVI formulations of the variably saturated water flow seepage problems in the fracture networks are proposed and a finite element procedure is also set up. To evaluate the unsaturated seepage flow in the heterogeneous fracture network, the following assumptions are made: the rock matrix is treated as impermeable; the fracture is nondeformable and the fracture system is under isothermal conditions; the fluid is essentially incompressible and it is assumed that the fluid flow obeys Darcy’s law. In contrast to the volume-averaged model, our seepage analysis is established on the basis of fracture segments and interconnections between fractures. The organization of the paper is as follows. In Section2, we define the governing equation of the partially saturated flow through single fractures with the form of modified Richards’ equation, and the complementary conditions of Signorini type on the seepage and infiltration surfaces are also presented. In Section 3, the proof of equivalence between the PDE and PVI formulations is provided. In Section 4, the corresponding numerical procedure is developed to solve the PVI formulations. In Section 5, four examples are discussed to demonstrate the reliability of the proposed algorithm and assess the sensitivity of variably saturated flow on the fracture distributions.


## 2. PDE Formulation for Variably Saturated Flow in Fracture Networks

### 2.1. Governing Equations

The variably saturated groundwater flow in a discrete fracture network, shown in Figure1, can be described by a governing equation for a single fracture. The fractures are usually conceptualized as parallel and smooth two-dimensional plates, in which the total head is uniform across the fracture width. The equation for variably saturated flow in a fracture of aperture bij can be derived by extending the saturated fracture flow equations and using the analogy of Richards’ equation [17] and the fact that the rate of change of water saturation is equal to the net inflow across the two endpoints of the fracture segment ij, including the net gain from fluid sources. With this extension, the governing equation for an arbitrary fracture segment ij,as shown in Figure 2, can be expressed in a modified form of Richards’ equation: (1)∂∂lbijkrkij∂ψ+z∂l=bijC∂ψ∂t,where ψ=ψl,t and z are the pressure head [L] and the elevation heads [L], respectively, in the fracture segment ij; kr is the relative permeability; C=∂Sij/∂ψ is the specific moisture capacity of fractures [L−1]; Sij is the water saturation. The saturated hydraulic conductivity kij [LT−1] of an individual fracture segment ij follows the cubic law, which can be stated as follows:(2)kij=gbij212υ,where g denotes the gravity acceleration [LT−2] and υ is the dynamic viscosity of water [L2T]; bij is the hydraulic fracture aperture [L], rather than the mean geometric aperture of the fracture.Figure 1
Variably saturated flow movement in the fracture network.Figure 2
Local coordinates of the fractureij.These two apertures are substantially different in general cases, except for smooth fractures. The cubic law, as the simplest approach to describe the fracture-dominated flow, was deduced on the basis of the fact that fracture was bounded by smooth and parallel plates. In reality, all fractures have rough walls and variable apertures, and fluid flow moving through a real fracture will take a tortuous path. Obviously, the deviation between real hydraulic properties and cubic law is expected [18, 19]. Brown [20] and Zimmerman and Bodvarsson [21] proposed mathematical expressions to relate these two apertures on the dependence of fracture roughness and contact area. When the fracture geometry information on fracture roughness and contact area is available or can be estimated, the fracture permeability can be evaluated through the cubic law as well as one of the relations between geometrical and hydraulic apertures. From the above discussion, an equivalent hydraulic aperture bij∗ [L] is commonly adopted to replace the mean aperture bij in (2).As shown in Figure3, intersections, the locations where two or more fractures are connected together, play an important role in conducting flow through a fracture network [22]. Supposing that there are mi fracture segments meeting at intersectioni,the sum of the total inflow flux and outflow flux at intersection i should be zero based on the principle of mass conservation, which can be expressed as(3)∑j=1mibijvij=0,where vij is an extension of Darcy flow velocity [LT−1] within the fracture segment ij under variable saturated flow condition, given by vij=-krkij∂ψ+z/∂l.Figure 3
Mass conservation at fracture intersectioni.


### 2.2. Constitutive Relationships

Constitutive relationships between the pressure headψ, relative permeability kr, and water saturation Sij are fundamental for modeling variably saturated underground flow processes. In order to solve the nonlinear flow equation (1) in terms of the pressure head, constitutive relationships must be available to describe variably saturated flow in the fractures. Unfortunately, there are a limited number of experimental, theoretical, and numerical studies, where these relationships were developed directly for variably saturated fracture flow.One approach to establish the related constitutive relationships for either single fractures [23] or fracture networks [4, 24] is to simply borrow or modify empirical functions from porous media. A commonly used functional relationship is the one proposed by van Genuchten [25], based on the earlier work by Mualem [26], where the pressure head-saturation and relative permeability-saturation relationships are expressed by (4)Se=1+αψn-mψ<01ψ≥0,kr=Se1/21-1-Se1/mm2,where α,n, and m are empirical fitting parameters with m=1-1/n, Se is the effective saturation given by Se=S-Sr/Ss-Sr, and Sr and Ss are the residual and satiated saturations, which are assigned to be 0 and 1 in this study, respectively.Nevertheless, the physical meaning of the relevant parameters in these relationships cannot always be easily understood for fractures, because the micro geometry in the fracture is two-dimensional while that of pores in porous media is three-dimensional [27]. The hydraulic characteristics of fractures are largely determined by fracture aperture and its connectivity. Currently, only a few studies have been conducted to directly quantify the variably saturated flow behavior as a function of fracture aperture [27–29]. Their results were based on capillary theory (percolation simulation) and they used the well-known cubic law to represent flow in the fracture. Recently, Ye et al. [30, 31] proposed simple and closed-form constitutive functions to describe the relationships between the relative permeability, saturation, and pressure head. The validity of these relationships was demonstrated by the consistency of laboratory test data with the numerical results of Chen and Horne [32] and Watanabe et al. (2015). On the basis of the aperture distribution function f(b), the relevant saturation and relative permeability yield(5)Se=∫bminbcbfbdb∫bmin∞bfbdbkr=Seβ∫bminbcb3fbdb∫bmin∞b3fbdb,where the local aperture larger than bc=2Tscos⁡α/ψ cannot be occupied by the water flow, Ts is the interfacial tension or surface tension, α is the contact angle, bmin is a critical aperture corresponding to the residual saturation, and β is an empirical parameter and takes the value of 1 in this study.


### 2.3. Boundary Conditions

Initial and boundary conditions are also needed to solve (1) describing unsaturated flow in the fractures. Boundary conditions can be on the pressure head, flux, seepage face, and infiltration boundaries.The seepage faceΓS=CD+DE shown in Figure 1 is a special external boundary of the saturated zone where water leaves the wet domain and ψ is uniformly zero, equal to the surrounding air pressure. When the seepage face is known a priori, it is considered as the pressure head boundary with a prescribed zero value, while the flux of the unsaturated zone qij is specified as zero with the negative pressure head. However, the positions of spill points between unsaturated and saturated zones are unknown a priori, and the singularity of spill points may result in mesh dependency and numerical instability [33] on (1). Based on the Signorini-type formulation, the potential seepage boundary condition ΓS for unsaturated seepage flow analysis in a fracture network can be described according to(6)qij≤0ψi=0i∈CDorqij=0ψi≤0i∈DE.On the rain infiltration surfacesΓR=EF+FG+GH+HI, the boundary condition is either a prescribed pressure head, as the rainwater infiltration rate is no more than the rainfall intensity, or a prescribed flux rate, as the pressure head is no more than the ponding depth. Similarly, the infiltration boundary can be defined as(7)qij≤bijIψi=ψpondi∈GH+HIorqij=bijIψi≤ψpondi∈EF+FG,where I is the rainfall intensity [LT−1] and ψpond is the ponding depth on the infiltration boundary [L].In nature, the conditions shown in (6) and (7) are also known as the unilateral boundary condition [34] (Hu et al., 2016) or the complementary condition of Signorini-type formulation [33, 35]. In the conditions, either the infiltration rate equals the intensity of rainfall or the water pressure head equals the ponding depth:(8)qij≤qu,ψi≤ψuqij-quψi-ψu=0i∈Γu=ΓS∪ΓR,where qu and ψu refer to critical flow rate [L2T−1] and water pressure head [L], respectively. From the above unified expression, we have qu=0,ψu=0 on the seepage boundary and qu=bijI,ψu=ψpond on the infiltration boundary.In addition, the saturated-unsaturated water flow in fracture networks should satisfy the following conditions:(1)
The initial condition:(9)ψx,z,0=ψ0x,z.(2)
The pressure head boundary condition(10)ψit=ψ-iti∈Γψ=AJ+BC.(3)
The flux boundary condition(11)qijt=bijvijt=q-iji∈Γq=AB,whereψ-i is the prescribed water head [L] and q-ij is the prescribed flow rate [L2T−1] for fracture ij.


## 2.1. Governing Equations

The variably saturated groundwater flow in a discrete fracture network, shown in Figure1, can be described by a governing equation for a single fracture. The fractures are usually conceptualized as parallel and smooth two-dimensional plates, in which the total head is uniform across the fracture width. The equation for variably saturated flow in a fracture of aperture bij can be derived by extending the saturated fracture flow equations and using the analogy of Richards’ equation [17] and the fact that the rate of change of water saturation is equal to the net inflow across the two endpoints of the fracture segment ij, including the net gain from fluid sources. With this extension, the governing equation for an arbitrary fracture segment ij,as shown in Figure 2, can be expressed in a modified form of Richards’ equation: (1)∂∂lbijkrkij∂ψ+z∂l=bijC∂ψ∂t,where ψ=ψl,t and z are the pressure head [L] and the elevation heads [L], respectively, in the fracture segment ij; kr is the relative permeability; C=∂Sij/∂ψ is the specific moisture capacity of fractures [L−1]; Sij is the water saturation. The saturated hydraulic conductivity kij [LT−1] of an individual fracture segment ij follows the cubic law, which can be stated as follows:(2)kij=gbij212υ,where g denotes the gravity acceleration [LT−2] and υ is the dynamic viscosity of water [L2T]; bij is the hydraulic fracture aperture [L], rather than the mean geometric aperture of the fracture.Figure 1
Variably saturated flow movement in the fracture network.Figure 2
Local coordinates of the fractureij.These two apertures are substantially different in general cases, except for smooth fractures. The cubic law, as the simplest approach to describe the fracture-dominated flow, was deduced on the basis of the fact that fracture was bounded by smooth and parallel plates. In reality, all fractures have rough walls and variable apertures, and fluid flow moving through a real fracture will take a tortuous path. Obviously, the deviation between real hydraulic properties and cubic law is expected [18, 19]. Brown [20] and Zimmerman and Bodvarsson [21] proposed mathematical expressions to relate these two apertures on the dependence of fracture roughness and contact area. When the fracture geometry information on fracture roughness and contact area is available or can be estimated, the fracture permeability can be evaluated through the cubic law as well as one of the relations between geometrical and hydraulic apertures. From the above discussion, an equivalent hydraulic aperture bij∗ [L] is commonly adopted to replace the mean aperture bij in (2).As shown in Figure3, intersections, the locations where two or more fractures are connected together, play an important role in conducting flow through a fracture network [22]. Supposing that there are mi fracture segments meeting at intersectioni,the sum of the total inflow flux and outflow flux at intersection i should be zero based on the principle of mass conservation, which can be expressed as(3)∑j=1mibijvij=0,where vij is an extension of Darcy flow velocity [LT−1] within the fracture segment ij under variable saturated flow condition, given by vij=-krkij∂ψ+z/∂l.Figure 3
Mass conservation at fracture intersectioni.


## 2.2. Constitutive Relationships

Constitutive relationships between the pressure headψ, relative permeability kr, and water saturation Sij are fundamental for modeling variably saturated underground flow processes. In order to solve the nonlinear flow equation (1) in terms of the pressure head, constitutive relationships must be available to describe variably saturated flow in the fractures. Unfortunately, there are a limited number of experimental, theoretical, and numerical studies, where these relationships were developed directly for variably saturated fracture flow.One approach to establish the related constitutive relationships for either single fractures [23] or fracture networks [4, 24] is to simply borrow or modify empirical functions from porous media. A commonly used functional relationship is the one proposed by van Genuchten [25], based on the earlier work by Mualem [26], where the pressure head-saturation and relative permeability-saturation relationships are expressed by (4)Se=1+αψn-mψ<01ψ≥0,kr=Se1/21-1-Se1/mm2,where α,n, and m are empirical fitting parameters with m=1-1/n, Se is the effective saturation given by Se=S-Sr/Ss-Sr, and Sr and Ss are the residual and satiated saturations, which are assigned to be 0 and 1 in this study, respectively.Nevertheless, the physical meaning of the relevant parameters in these relationships cannot always be easily understood for fractures, because the micro geometry in the fracture is two-dimensional while that of pores in porous media is three-dimensional [27]. The hydraulic characteristics of fractures are largely determined by fracture aperture and its connectivity. Currently, only a few studies have been conducted to directly quantify the variably saturated flow behavior as a function of fracture aperture [27–29]. Their results were based on capillary theory (percolation simulation) and they used the well-known cubic law to represent flow in the fracture. Recently, Ye et al. [30, 31] proposed simple and closed-form constitutive functions to describe the relationships between the relative permeability, saturation, and pressure head. The validity of these relationships was demonstrated by the consistency of laboratory test data with the numerical results of Chen and Horne [32] and Watanabe et al. (2015). On the basis of the aperture distribution function f(b), the relevant saturation and relative permeability yield(5)Se=∫bminbcbfbdb∫bmin∞bfbdbkr=Seβ∫bminbcb3fbdb∫bmin∞b3fbdb,where the local aperture larger than bc=2Tscos⁡α/ψ cannot be occupied by the water flow, Ts is the interfacial tension or surface tension, α is the contact angle, bmin is a critical aperture corresponding to the residual saturation, and β is an empirical parameter and takes the value of 1 in this study.


## 2.3. Boundary Conditions

Initial and boundary conditions are also needed to solve (1) describing unsaturated flow in the fractures. Boundary conditions can be on the pressure head, flux, seepage face, and infiltration boundaries.The seepage faceΓS=CD+DE shown in Figure 1 is a special external boundary of the saturated zone where water leaves the wet domain and ψ is uniformly zero, equal to the surrounding air pressure. When the seepage face is known a priori, it is considered as the pressure head boundary with a prescribed zero value, while the flux of the unsaturated zone qij is specified as zero with the negative pressure head. However, the positions of spill points between unsaturated and saturated zones are unknown a priori, and the singularity of spill points may result in mesh dependency and numerical instability [33] on (1). Based on the Signorini-type formulation, the potential seepage boundary condition ΓS for unsaturated seepage flow analysis in a fracture network can be described according to(6)qij≤0ψi=0i∈CDorqij=0ψi≤0i∈DE.On the rain infiltration surfacesΓR=EF+FG+GH+HI, the boundary condition is either a prescribed pressure head, as the rainwater infiltration rate is no more than the rainfall intensity, or a prescribed flux rate, as the pressure head is no more than the ponding depth. Similarly, the infiltration boundary can be defined as(7)qij≤bijIψi=ψpondi∈GH+HIorqij=bijIψi≤ψpondi∈EF+FG,where I is the rainfall intensity [LT−1] and ψpond is the ponding depth on the infiltration boundary [L].In nature, the conditions shown in (6) and (7) are also known as the unilateral boundary condition [34] (Hu et al., 2016) or the complementary condition of Signorini-type formulation [33, 35]. In the conditions, either the infiltration rate equals the intensity of rainfall or the water pressure head equals the ponding depth:(8)qij≤qu,ψi≤ψuqij-quψi-ψu=0i∈Γu=ΓS∪ΓR,where qu and ψu refer to critical flow rate [L2T−1] and water pressure head [L], respectively. From the above unified expression, we have qu=0,ψu=0 on the seepage boundary and qu=bijI,ψu=ψpond on the infiltration boundary.In addition, the saturated-unsaturated water flow in fracture networks should satisfy the following conditions:(1)
The initial condition:(9)ψx,z,0=ψ0x,z.(2)
The pressure head boundary condition(10)ψit=ψ-iti∈Γψ=AJ+BC.(3)
The flux boundary condition(11)qijt=bijvijt=q-iji∈Γq=AB,whereψ-i is the prescribed water head [L] and q-ij is the prescribed flow rate [L2T−1] for fracture ij.


## 3. PVI Formulation for Variably Saturated Flow in Fracture Networks

Richards’ equation for describing variably saturated flow in the fracture networks is nonlinear and the solution is further complicated due to the singularity of seepage face and rainfall boundary conditions, as shown in (8). Thus, it is essential to develop an efficient and robust approach to obtain an accurate solution to the governing equations. The PVI approach based on a rigorous mathematical foundation can weaken the nonlinearity caused by the unilateral boundary condition, in which the flux part is transformed into the natural boundary condition. In this study, the PVI formulation of the variably saturated flow in fracture networks is stated below: for a specified set of trial functions(12)ΦPVI=ψ∣ψi=ψ-onΓψ;ψi≤ψuonΓu,the theoretical solution ψ∈ΦPVI with respect to ∀φ∈ΦPVI should satisfy(13)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ≥0,where (14)∂ψ∂t,φ-ψ=∑Ω∫lijφ-ψbijC∂ψ∂tdl(15)ξψ,φ-ψ=∑Ω∫lij∂φ-ψ∂lbijkrkij∂ψ+z∂ldl(16)ζφ-ψ=∑i∈Γq∑j=1mi-φi-ψiq-ij+∑i∈Γu∑j=1mi-φi-ψiqu.By employing integration by parts, (15) is given by(17)ξψ,φ-ψ=∑Ω∫lij∂φ-ψ∂lbijkrkij∂ψ+z∂ldl=∑Ω∫lijφ-ψdl-φ-ψbijvij0lij=∑Ω∫lijφ-ψ∂∂lbijkrkij∂ψ+z∂ldl+∑i=1n∑j=1miφi-ψibijvij,where n is the total number of fracture nodes.Inserting (13), (16), and (17) into (11), the following relationship can be obtained:(18)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑Ω∫lijφ-ψ∂∂lbijkrkij∂ψ+z∂l+bijC∂ψ∂tdl+∑i∈Ωin∑j=1miφi-ψibijvij+∑i∈Γq∑j=1miφi-ψiqij-q-ij+∑i∈Γu∑j=1miφi-ψiqij-qu⩾0.


### 3.1. The Proof for PDE⇒ PVI

Suppose that the function ofψ is the solution for the PDE formulations. Combining (1), (3), and (11) and the fact that ψ belongs to ΦPVI gives(19)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Γu∑j=1miφi-ψiqij-qu=∑i∈Γu∑j=1miφi-ψuqij-qu-∑i∈Γu∑j=1miψi-ψuqij-qu=∑i∈Γu∑j=1miφi-ψuqij-qu≥0.Hence,ψ is the solution for the PVI.


### 3.2. The Proof for PVI⇒ PDE

Similarly,ψ is assumed as the solution for the PVI, and (18) holds for ∀φ∈ΦPVI. Inserting φ=ψ+ε1 and φ=ψ-ε1 into (18), where ε1 is an arbitrary function whose values are zeros at all the nodes in the whole fracture network, the continuity equation (1) can be derived:(20)∂∂lbijkrkij∂ψ+z∂l=bijC∂ψ∂t.Thus, (18) can be expressed as(21)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Ωin∑j=1miφi-ψibijvij+∑i∈Γq∑j=1miφi-ψiqij-q-ij+∑i∈Γu∑j=1miφi-ψiqij-qu⩾0.Supposingφ=ψ+ε2 and φ=ψ-ε2 in (19), where ε2 is any function that becomes zero at the nodes on Γψ, Γq, and Γu, the mass conservation equation on the intersections can be obtained:(22)∑j=1mibijvij=0.Equation (21) can be rewritten as(23)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Γq∑j=1miφi-ψiqij-q-ij+∑i∈Γu∑j=1miφi-ψiqij-qu⩾0.Takingφ=ψ+ε3 and φ=ψ-ε3 in (21), with ε3 being any function that equals zero at the nodes of boundaries Γψ and Γu, the flux condition yields(24)qijt=q-ij.Herein, (23) can be given as(25)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Γu∑j=1miφi-ψuqij-qu-∑i∈Γu∑j=1miψi-ψuqij-qu⩾0.As the term∑i∈Γu∑j=1miψi-ψuqij-qu in the above inequality is a definite value, the inequality condition of (23) for ∀φ∈ΦPVI leads to(26)∑i∈Γu∑j=1miφi-ψuqij-qu⩾0.From (12), φ is an arbitrary function that satisfies φi≤ψu for Γu, so the following equation is obtained:(27)qij-qu≤0.Particularly, insertingφ=ψu in (23) results in(28)∑i∈Γu∑j=1miψi-ψuqij-qu≤0.From (27) and (28) and the fact that ψi≤ψu for Γu, the complementary condition of Signorini-type formulation for Γu is expressed as(29)qij-quψi-ψu=0.The equivalence relation between the PDE and PVI formulations is presented in the above proof. It is important to note that the PVI formulation as shown in (13) can naturally satisfy the complementary condition on Γu, and the difficulty in selecting the trial functions arising from the boundary nonlinearity is also reduced during the solving process.


## 3.1. The Proof for PDE⇒ PVI

Suppose that the function ofψ is the solution for the PDE formulations. Combining (1), (3), and (11) and the fact that ψ belongs to ΦPVI gives(19)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Γu∑j=1miφi-ψiqij-qu=∑i∈Γu∑j=1miφi-ψuqij-qu-∑i∈Γu∑j=1miψi-ψuqij-qu=∑i∈Γu∑j=1miφi-ψuqij-qu≥0.Hence,ψ is the solution for the PVI.


## 3.2. The Proof for PVI⇒ PDE

Similarly,ψ is assumed as the solution for the PVI, and (18) holds for ∀φ∈ΦPVI. Inserting φ=ψ+ε1 and φ=ψ-ε1 into (18), where ε1 is an arbitrary function whose values are zeros at all the nodes in the whole fracture network, the continuity equation (1) can be derived:(20)∂∂lbijkrkij∂ψ+z∂l=bijC∂ψ∂t.Thus, (18) can be expressed as(21)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Ωin∑j=1miφi-ψibijvij+∑i∈Γq∑j=1miφi-ψiqij-q-ij+∑i∈Γu∑j=1miφi-ψiqij-qu⩾0.Supposingφ=ψ+ε2 and φ=ψ-ε2 in (19), where ε2 is any function that becomes zero at the nodes on Γψ, Γq, and Γu, the mass conservation equation on the intersections can be obtained:(22)∑j=1mibijvij=0.Equation (21) can be rewritten as(23)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Γq∑j=1miφi-ψiqij-q-ij+∑i∈Γu∑j=1miφi-ψiqij-qu⩾0.Takingφ=ψ+ε3 and φ=ψ-ε3 in (21), with ε3 being any function that equals zero at the nodes of boundaries Γψ and Γu, the flux condition yields(24)qijt=q-ij.Herein, (23) can be given as(25)∂ψ∂t,φ-ψ+ξψ,φ-ψ+ζφ-ψ=∑i∈Γu∑j=1miφi-ψuqij-qu-∑i∈Γu∑j=1miψi-ψuqij-qu⩾0.As the term∑i∈Γu∑j=1miψi-ψuqij-qu in the above inequality is a definite value, the inequality condition of (23) for ∀φ∈ΦPVI leads to(26)∑i∈Γu∑j=1miφi-ψuqij-qu⩾0.From (12), φ is an arbitrary function that satisfies φi≤ψu for Γu, so the following equation is obtained:(27)qij-qu≤0.Particularly, insertingφ=ψu in (23) results in(28)∑i∈Γu∑j=1miψi-ψuqij-qu≤0.From (27) and (28) and the fact that ψi≤ψu for Γu, the complementary condition of Signorini-type formulation for Γu is expressed as(29)qij-quψi-ψu=0.The equivalence relation between the PDE and PVI formulations is presented in the above proof. It is important to note that the PVI formulation as shown in (13) can naturally satisfy the complementary condition on Γu, and the difficulty in selecting the trial functions arising from the boundary nonlinearity is also reduced during the solving process.


## 4. Numerical Solution

When the finite element method and a backward time difference scheme are applied to (13), each fracture segment is treated as a line element in the whole fracture network domain. The finite element discretization of the PVI formulation (13) for variably saturated water flow in a discrete fracture network can be presented as follows: find a vector ψη+1∈ΦPVIψ corresponding to ∀φ∈ΦPVIψ in a trial vector space(30)ΦPVIψ=ψ∣ψ∈Rn;ψi=ψ-ifori∈Γψ;ψi≤ψufori∈Γuwhich should obey the following relationships:(31)φ-ψη+1TC+ΔtKψη+1≥φ-ψη+1TCψη+φ-ψη+1TΔtq,where(32)K=∑Ω∫lijBTbijkrkijBdlC=∑Ω∫lijNTbijCNdlq=∑i∈ΓqNTqij+∑i∈ΓuNTqu-∑Ω∫lijBTbijkrkijBdl·zN=1-llijllij,B=-1lij1lij,where η is the time step, Δt is the time increment, and N and B are the matrices of shape and geometric functions.Equation (31) is the two-dimensional finite element algorithm for approximating the PVI formulation by matrix representation. The numerical procedure for solving (31) involves a time-marching loop and a nested loop to confirm the seepage face and infiltration boundaries. For the time-marching loop, a constant time increment is adopted with a small size to guarantee the accuracy of numerical simulation. For the nested loop at each time step, the complementary algorithm proposed by Jiang et al. [15, 36] is employed to eliminate the singularity of the seepage face and infiltration boundaries. The convergence criterion for the finite element procedure in this study is given as follows:(33)ψτη-ψτ-1η⩽δψτ-1η-ψη-1,where τ is the iteration step and the symbol δ denotes user-specified error tolerance and takes the value of 0.001 in this study.


## 5. Verification and Illustrative Examples

Four typical examples are presented to investigate the performance of the proposed algorithm. In the first and second examples, the van Genuchten [25] model developed from Mualem’s model [26] is employed during numerical simulation; the relevant parameters are Ss=0.363,Sr=0.186,α=1.0 m−1, n=1.53, and kij=1.0×10-6orbij=1.11×10-6. The other two examples are based on the work of Ye et al. [16, 31] and the relevant parameters are (1)  bmin=0,β=1.5 and (2)  bmin=0,β=1. The three related retention curves are presented in Figure 4.Figure 4
Constitutive relationships for four examples in Section5: (a) pressure head versus saturation relationships and (b) relative permeability versus saturation relationships.
(a)
(b)


### 5.1. One-Dimensional Infiltration Problems

The first example validates the proposed method with respect to the semianalytical solution derived by Warrick et al. [37], where the one-dimensional infiltration problem in homogeneous media was considered. As shown in Figure 5, a vertical fracture with a constant aperture of 1.11 mm is modeled with 1 m height, and an initial pressure head of −8 m is assigned to the whole fracture domain. A total of 50 line elements with a size of 0.02 m and 51 nodes are generated to mesh the domain, and the time step size keeps constant at 0.05 hr. To define the boundary conditions, a zero pressure head is specified at the fracture top for seepage analysis, while the bottom of the fracture still keeps -8 m pressure head.Figure 5
Simulated and analytical water table positions at various times.According to the transient analysis through the one-dimensional infiltration fracture, comparisons of the numerical predictions and Warrick’s solutions at three elapsed times are shown in Figure5. The calculated results agree well with the theoretical data, indicating that the proposed method on the basis of PVI formulations can reasonably represent the variably saturated flow through the fractures. Readers may refer to Warrick et al. [37] for more details about the analytical solution for the one-dimensional infiltration problems in homogeneous media.


### 5.2. Steady Flow in Homogeneous and Nonhomogenous Rectangular Dams

Example2 involves unconfined seepage through homogeneous and nonhomogenous rectangular dams driven by the head difference between 10 m at the upstream face and 2 m at the downstream face (Figure 6). The absolute permeability k of the homogenous rectangular dam (Figure 6(a)) is 1 m/d. The nonhomogenous rectangular dam contains two different materials with a material boundary in the middle (Figure 6(b)). The absolute permeability to the right of the material boundary is 10 times higher than that to the left. We use our developed PVI model to simulate this seepage problem and locate the steady-state free surface. Two equivalent fracture networks formed by two orthogonal sets of fractures are used to predict the free surface of homogeneous and nonhomogenous rectangular dams. The aperture values corresponding to 1 m/d and 10 m/d are equal to 0.121 mm and 0.261 mm according to the principle of flux equivalence [16, 31].Figure 6
Verification on saturated steady seepage problems.
(a)
Homogenous dam
(b)
Nonhomogenous damFigure6 shows the results for the steady-state free surfaces predicted by the proposed PVI method. In comparison, the numerical results performed by Lacy and Prevost [38] and Borja and Kishnani [39] are also presented in Figure 6. It can be seen that the results of the free surfaces in the nonhomogenous dam have good agreements, but the results of the free surfaces in the homogenous dam calculated by the proposed PVI method are closer to the numerical data obtained by Lacy and Prevost [38] than those provided by Borja and Kishnani [39]. Note that the seepage point calculated by Borja and Kishnani [39] has a very pronounced difference compared to the proposed PVI method and Lacy and Prevost [38], which results from the singularity of the seepage face. In the proposed PVI method, the singularity of the seepage face is transformed into the natural boundary condition to guarantee the numerical stability and validity.


### 5.3. Preferential Flow Model

The third example involves groundwater infiltrations in two typical 1.0 m × 1.0 m squares with simple two-dimensional fracture networks consisting of vertical and horizontal fractures, as shown in Figure7. We assume the nonflow conditions at the right and left side boundaries. The water pressure head at the top boundary remains constant at its initial value of 1 m, whereas the bottom boundary is imposed as Signorini-type complementary condition. The initial water pressure head is specified as −1 m.Figure 7
Numerical generations of two fracture systems: (a) fracture sets with mean trace length of 2 m and spacing of 0.005 m and (b) fracture sets with mean trace length of 0.2 m and spacing of 0.01 m.
(a)
Case1
(b)
Case2Different methods exist to computationally generate a fracture network [40], and the discrete fracture networks in this study are constructed by the Monte-Carlo method [41]. Note that all or most connected fractures are expected to participate in variably saturated flow in the fracture network and that the fracture dead ends and fractures that are isolated from the major part of the network may not be hydraulically connected and thus can be disregarded in studies of unsaturated flow processes in fracture networks. Similar fracture networks have also been implemented by other researchers to study flow and transport properties of fractures [24, 42]To investigate how the extent and density of fractures affect the variably saturated flow infiltration processes in the rock square, two types of fracture systems are mimicked and their related statistical parameters are given in Table1. Obviously, the horizontal and vertical fractures have the same statistical parameters in each case, except the aperture size. However, the mean trace length in case 1 is ten times that in case 2, and the mean spacing in case 1 is half that in case 2. The corresponding realizations of two different fracture systems are presented in Figure 7. For numerical modeling, a time step of 1 s is adopted.Table 1
Parameters of fractures and probability model.
Dip (°)
Aperture (10−4 m)
Mean trace length (m)
Mean spacing (m)
Mean values
Variance
Probability model
Case1
180
0.1
0.0
Normal
2.0
0.005
90
0.5
0.5
Normal
2.0
0.005
Case2
180
0.1
0.0
Normal
0.2
0.010
90
0.5
0.5
Normal
0.2
0.010Figure8 shows the simulated seeps in the two heterogeneous fracture systems for the top boundary pressure head of 1 m. As water is injected at a prescribed pressure head into the fractures, water pressure head increases and inflow into the fracture network beneath the top boundary takes place. The inflow spreads nonuniformly through the heterogeneous fractures on the basis of spatially variable aperture distribution and unsaturated condition.Figure 8
Variation of water seepage in the two fracture systems at times 2 s and 10 s.
(a)
Case1
(b)
Case2The seepage patterns are distinguishable in appearance but share some common features. Flow generally proceeds in a manner of narrow fingers in the early stage. Several flow paths can develop from localized infiltration. As indicated by Pruess [12], when the fast preferential flow is driven by gravity force in the downward direction in a continuum model, the seepage velocity can be represented as a function of the volumetric flux, porosity, and water saturation. The fast preferential flow can occur in regions with larger Darcy velocity, smaller porosity, and smaller effective water saturation. According to the cubic law, the fractures with larger aperture values can give rise to higher localized Darcy velocity resulting from the heterogeneities of the fracture distribution. Hence, influenced by the random distribution of fracture aperture, the fast flow pathways are dominated by the factures with larger aperture values as displayed in Figure 8.Comparing cases1 and 2, it is obvious that the vertical advancements of water infiltration are substantially different. Fractures in Figure 8(a) have a greater trace length and density, and the evolution of water seeps can develop faster along vertical fractures and the breakthrough on the bottom boundary can take place in a shorter time. As a result, the downward advancement of seepage flow in Figure 8(b) is fairly slower than in Figure 8(a), due to the connectivity of fracture networks.The connectivity of the fracture networks is an essential feature controlling the flow movement in these impermeable geological fractured media. In order to quantify the connectivity in these heterogeneous fracture networks with fracture geometrical properties, the concept of geological entropy [43] is employed to characterize the heterogeneity of these fracture networks rather than the fracture density (the number of fractures per unit area). Based on the theoretical and numerical study of the geological entropy in porous media, the related relationships yield(34)HR=H-LHG(35)HG=-∑n=1NpG,nln⁡pG,n(36)H-L=1I×J∑j=1J∑i=1IHLi,j(37)HLl=-∑n=1Npl,nln⁡pl,n,where HR is the relative entropy index, which is a quantitative measure of the spatial disorder in a discrete random system; HG is the global entropy of the discrete random system representing a set of N grid blocks; HL is the local entropy of each square grid; pG’s are defined as the volumetric proportions of the fractures or rock masses over the whole domain of interest; pl’s are defined as the volumetric proportions of the fractures or rock masses over each square grid; H-L is the average of the local entropy over the entire domain.An example to characterize such fracture network system is shown in Figure9; two groups of 1.0 m × 1.0 m squares are composed of two simple vertical and horizontal fractures. The lengths of apertures for both systems are 0.5 m and 0.01 m, respectively, but the positions are different. It can be observed intuitively that the left plot domain is unconnected or disordered while the right is connected or ordered. In both systems, pG’s of fractures and rock masses are equal to 0.02 and 0.98 according to(35), corresponding to an HG value of 0.098. When the two systems are divided into nine subsets as shown in Figure 9, the relative entropy HR values are reported in Table 2. The HR index for the connected system is equal to 0.867. On the other hand, the unconnected system has a very large HR value, which is 0.969. It is found that the indicator HR is sensitive to the connectivity of the fracture networks. Furthermore, Bianchi and Pedretti [43] pointed out that HR is a global metric that integrates multiple properties of the spatial distribution of hydraulic conductivity.Table 2
Parameters of the geological entropy for connected and unconnected systems.
Subdomain
1
2
3
4
5
6
7
8
9
Connected system
pl
Fracture
0
0.03
0
0.03
0.06
0.03
0
0.03
0
Rock mass
1
0.97
1
0.97
0.94
0.97
1
0.97
1
Hl
0
0.135
0
0.135
0.227
0.135
0
0.135
0
H-L
0.085
HG
0.098
HR
0.867
Unconnected system
pl
Fracture
0.03
0.015
0.03
0.015
0
0.015
0.03
0.015
0.03
Rock mass
0.97
0.985
0.97
0.985
1
0.985
0.97
0.985
0.97
Hl
0.135
0.078
0.135
0.078
0
0.078
0.135
0.078
0.135
H-L
0.095
HG
0.098
HR
0.969Figure 9
Local and relative entropy calculation for connected and unconnected systems.
(a)
Connected (HR=0.876)
(b)
Unconnected (HR=0.969)Herein, the relative entropyHR is also imposed to measure the connectivity of the fracture distribution plotted in Figure 7 with different sizes of square grid, which is shown in Figure 10. In both models, HR increases with the size of square grid. However, the HR value in case 1 with faster finger flow is much smaller than that in case 2 with a slower volume-averaged manner. The results are consistent with the conclusion of Bianchi and Pedretti [43] that fingering is more pronounced when the relative entropy HR is low. Therefore, in a way similar to the derivation of geological entropy indicators, it can be recognized that the relative entropy HR, which takes into account the spatial distributions of the heterogeneous fracture network including connectivity, extent and spacing, and aperture distribution, has significant effects on the development and propagation of preferential flow during heterogeneous seepage analysis. Note that the relative entropy HR cannot be simply treated as an alternative to connectivity, but rather a comprehensive description including the connectivity, the spatial distribution of aperture, and extent and orientation.Figure 10
Variation ofHR with different square grid size.


### 5.4. Influences of Fracture Orientation on Rainfall Infiltration

The fourth example is to investigate the effect of rainfall infiltration on a fractured rock slope with two sets of random distribution fractures, where the rapid generation of higher hydraulic pressure within a short possible time may become one of the critical factors on slope stability. Four types of discrete fracture networks are employed and their geometrical parameters are given in Table3. The corresponding computational meshes are created using the Monte-Carlo method as presented in Figure 11. Suppose that the strong rainfall process lasts for 50 hours with a constant intensity of 0.001 m/hr. In the analysis, the time increment of each time step is adopted as 1.0 hr.Table 3
Parameters of fractures and probability model.
Dip (°)
Aperture (10−4 m)
Mean trace length (m)
Mean spacing (m)
Mean values
Variance
Probability model
Case1
0.1
Normal
0.2
15
0.3
2.0
135
0.1
0.8
Case2
30
0.3
2.0
135
0.1
0.8
Case3
45
0.3
2.0
135
0.1
0.8
Case4
60
0.3
2.0
135
0.1
0.8Figure 11
Numerical realizations of two fracture systems in a fractured rock slope under rainfall condition: (a) dip of 15 degrees; (b) dip of 30 degrees; (c) dip of 45 degrees; (d) dip of 60 degrees.
(a)
Case1
(b)
Case2
(c)
Case3
(d)
Case4The left and bottom boundaries are impermeable and the right boundary under the initial water table is always treated as a water pressure head boundary. The infiltration boundaries are imposed on the other boundaries along the slope. Due to the high gradient of the slope, it is reasonable to assume that no runoff would be produced on the infiltration boundary, which creates a zero value ofψpond. As an initial condition, we assume that the domain is under a hydrostatic state with an external water table situated at 5.00 m. The water pressure head below the groundwater level increases linearly from 0 m at the water table position to 5 m at the bottom, and that above the groundwater level decreases linearly from 0 m at the water table position to -20 m at the top of the slope.The normalized flow rate distributions inside the rock slope after 10 and 20 hours of water infiltration are shown in Figure12. The seepage flow occurs prevalently along bedding plane fractures caused by higher hydraulic conductivity of larger facture aperture; a portion of factures are still not active during rainfall infiltration. It can be observed that, with the increase of fracture dip, the propagation of water seepage extends inward more rapidly and the expansion of transient saturated region is driven faster into a deeper domain. The fractures with a steep dip can give rise to fast movement of water infiltration. A comparison of the normalized flow rate distributions ranging from 10 to 20 hr indicates that the orientation of fractures is also an important feature for describing the evolution of fast seepage flow. Figure 13 shows the variation of the total net flux into the unsaturated domain with time. When the fracture dip increases from 15 to 45 degrees, the growth of the total net flux increases substantially; conversely, the growth of the total net flux turns out to be smaller when the fracture dip reaches 60 degrees, especially after 20 hours water infiltration.Figure 12
The normalized flow rate distributions in the slope at times (a) 10 hr and (b) 20 hr.
(a)
Case1
(b)
Case2
(c)
Case3
(d)
Case4Figure 13
Evolution of the total net flow into the unsaturated domain.


## 5.1. One-Dimensional Infiltration Problems

The first example validates the proposed method with respect to the semianalytical solution derived by Warrick et al. [37], where the one-dimensional infiltration problem in homogeneous media was considered. As shown in Figure 5, a vertical fracture with a constant aperture of 1.11 mm is modeled with 1 m height, and an initial pressure head of −8 m is assigned to the whole fracture domain. A total of 50 line elements with a size of 0.02 m and 51 nodes are generated to mesh the domain, and the time step size keeps constant at 0.05 hr. To define the boundary conditions, a zero pressure head is specified at the fracture top for seepage analysis, while the bottom of the fracture still keeps -8 m pressure head.Figure 5
Simulated and analytical water table positions at various times.According to the transient analysis through the one-dimensional infiltration fracture, comparisons of the numerical predictions and Warrick’s solutions at three elapsed times are shown in Figure5. The calculated results agree well with the theoretical data, indicating that the proposed method on the basis of PVI formulations can reasonably represent the variably saturated flow through the fractures. Readers may refer to Warrick et al. [37] for more details about the analytical solution for the one-dimensional infiltration problems in homogeneous media.


## 5.2. Steady Flow in Homogeneous and Nonhomogenous Rectangular Dams

Example2 involves unconfined seepage through homogeneous and nonhomogenous rectangular dams driven by the head difference between 10 m at the upstream face and 2 m at the downstream face (Figure 6). The absolute permeability k of the homogenous rectangular dam (Figure 6(a)) is 1 m/d. The nonhomogenous rectangular dam contains two different materials with a material boundary in the middle (Figure 6(b)). The absolute permeability to the right of the material boundary is 10 times higher than that to the left. We use our developed PVI model to simulate this seepage problem and locate the steady-state free surface. Two equivalent fracture networks formed by two orthogonal sets of fractures are used to predict the free surface of homogeneous and nonhomogenous rectangular dams. The aperture values corresponding to 1 m/d and 10 m/d are equal to 0.121 mm and 0.261 mm according to the principle of flux equivalence [16, 31].Figure 6
Verification on saturated steady seepage problems.
(a)
Homogenous dam
(b)
Nonhomogenous damFigure6 shows the results for the steady-state free surfaces predicted by the proposed PVI method. In comparison, the numerical results performed by Lacy and Prevost [38] and Borja and Kishnani [39] are also presented in Figure 6. It can be seen that the results of the free surfaces in the nonhomogenous dam have good agreements, but the results of the free surfaces in the homogenous dam calculated by the proposed PVI method are closer to the numerical data obtained by Lacy and Prevost [38] than those provided by Borja and Kishnani [39]. Note that the seepage point calculated by Borja and Kishnani [39] has a very pronounced difference compared to the proposed PVI method and Lacy and Prevost [38], which results from the singularity of the seepage face. In the proposed PVI method, the singularity of the seepage face is transformed into the natural boundary condition to guarantee the numerical stability and validity.


## 5.3. Preferential Flow Model

The third example involves groundwater infiltrations in two typical 1.0 m × 1.0 m squares with simple two-dimensional fracture networks consisting of vertical and horizontal fractures, as shown in Figure7. We assume the nonflow conditions at the right and left side boundaries. The water pressure head at the top boundary remains constant at its initial value of 1 m, whereas the bottom boundary is imposed as Signorini-type complementary condition. The initial water pressure head is specified as −1 m.Figure 7
Numerical generations of two fracture systems: (a) fracture sets with mean trace length of 2 m and spacing of 0.005 m and (b) fracture sets with mean trace length of 0.2 m and spacing of 0.01 m.
(a)
Case1
(b)
Case2Different methods exist to computationally generate a fracture network [40], and the discrete fracture networks in this study are constructed by the Monte-Carlo method [41]. Note that all or most connected fractures are expected to participate in variably saturated flow in the fracture network and that the fracture dead ends and fractures that are isolated from the major part of the network may not be hydraulically connected and thus can be disregarded in studies of unsaturated flow processes in fracture networks. Similar fracture networks have also been implemented by other researchers to study flow and transport properties of fractures [24, 42]To investigate how the extent and density of fractures affect the variably saturated flow infiltration processes in the rock square, two types of fracture systems are mimicked and their related statistical parameters are given in Table1. Obviously, the horizontal and vertical fractures have the same statistical parameters in each case, except the aperture size. However, the mean trace length in case 1 is ten times that in case 2, and the mean spacing in case 1 is half that in case 2. The corresponding realizations of two different fracture systems are presented in Figure 7. For numerical modeling, a time step of 1 s is adopted.Table 1
Parameters of fractures and probability model.
Dip (°)
Aperture (10−4 m)
Mean trace length (m)
Mean spacing (m)
Mean values
Variance
Probability model
Case1
180
0.1
0.0
Normal
2.0
0.005
90
0.5
0.5
Normal
2.0
0.005
Case2
180
0.1
0.0
Normal
0.2
0.010
90
0.5
0.5
Normal
0.2
0.010Figure8 shows the simulated seeps in the two heterogeneous fracture systems for the top boundary pressure head of 1 m. As water is injected at a prescribed pressure head into the fractures, water pressure head increases and inflow into the fracture network beneath the top boundary takes place. The inflow spreads nonuniformly through the heterogeneous fractures on the basis of spatially variable aperture distribution and unsaturated condition.Figure 8
Variation of water seepage in the two fracture systems at times 2 s and 10 s.
(a)
Case1
(b)
Case2The seepage patterns are distinguishable in appearance but share some common features. Flow generally proceeds in a manner of narrow fingers in the early stage. Several flow paths can develop from localized infiltration. As indicated by Pruess [12], when the fast preferential flow is driven by gravity force in the downward direction in a continuum model, the seepage velocity can be represented as a function of the volumetric flux, porosity, and water saturation. The fast preferential flow can occur in regions with larger Darcy velocity, smaller porosity, and smaller effective water saturation. According to the cubic law, the fractures with larger aperture values can give rise to higher localized Darcy velocity resulting from the heterogeneities of the fracture distribution. Hence, influenced by the random distribution of fracture aperture, the fast flow pathways are dominated by the factures with larger aperture values as displayed in Figure 8.Comparing cases1 and 2, it is obvious that the vertical advancements of water infiltration are substantially different. Fractures in Figure 8(a) have a greater trace length and density, and the evolution of water seeps can develop faster along vertical fractures and the breakthrough on the bottom boundary can take place in a shorter time. As a result, the downward advancement of seepage flow in Figure 8(b) is fairly slower than in Figure 8(a), due to the connectivity of fracture networks.The connectivity of the fracture networks is an essential feature controlling the flow movement in these impermeable geological fractured media. In order to quantify the connectivity in these heterogeneous fracture networks with fracture geometrical properties, the concept of geological entropy [43] is employed to characterize the heterogeneity of these fracture networks rather than the fracture density (the number of fractures per unit area). Based on the theoretical and numerical study of the geological entropy in porous media, the related relationships yield(34)HR=H-LHG(35)HG=-∑n=1NpG,nln⁡pG,n(36)H-L=1I×J∑j=1J∑i=1IHLi,j(37)HLl=-∑n=1Npl,nln⁡pl,n,where HR is the relative entropy index, which is a quantitative measure of the spatial disorder in a discrete random system; HG is the global entropy of the discrete random system representing a set of N grid blocks; HL is the local entropy of each square grid; pG’s are defined as the volumetric proportions of the fractures or rock masses over the whole domain of interest; pl’s are defined as the volumetric proportions of the fractures or rock masses over each square grid; H-L is the average of the local entropy over the entire domain.An example to characterize such fracture network system is shown in Figure9; two groups of 1.0 m × 1.0 m squares are composed of two simple vertical and horizontal fractures. The lengths of apertures for both systems are 0.5 m and 0.01 m, respectively, but the positions are different. It can be observed intuitively that the left plot domain is unconnected or disordered while the right is connected or ordered. In both systems, pG’s of fractures and rock masses are equal to 0.02 and 0.98 according to(35), corresponding to an HG value of 0.098. When the two systems are divided into nine subsets as shown in Figure 9, the relative entropy HR values are reported in Table 2. The HR index for the connected system is equal to 0.867. On the other hand, the unconnected system has a very large HR value, which is 0.969. It is found that the indicator HR is sensitive to the connectivity of the fracture networks. Furthermore, Bianchi and Pedretti [43] pointed out that HR is a global metric that integrates multiple properties of the spatial distribution of hydraulic conductivity.Table 2
Parameters of the geological entropy for connected and unconnected systems.
Subdomain
1
2
3
4
5
6
7
8
9
Connected system
pl
Fracture
0
0.03
0
0.03
0.06
0.03
0
0.03
0
Rock mass
1
0.97
1
0.97
0.94
0.97
1
0.97
1
Hl
0
0.135
0
0.135
0.227
0.135
0
0.135
0
H-L
0.085
HG
0.098
HR
0.867
Unconnected system
pl
Fracture
0.03
0.015
0.03
0.015
0
0.015
0.03
0.015
0.03
Rock mass
0.97
0.985
0.97
0.985
1
0.985
0.97
0.985
0.97
Hl
0.135
0.078
0.135
0.078
0
0.078
0.135
0.078
0.135
H-L
0.095
HG
0.098
HR
0.969Figure 9
Local and relative entropy calculation for connected and unconnected systems.
(a)
Connected (HR=0.876)
(b)
Unconnected (HR=0.969)Herein, the relative entropyHR is also imposed to measure the connectivity of the fracture distribution plotted in Figure 7 with different sizes of square grid, which is shown in Figure 10. In both models, HR increases with the size of square grid. However, the HR value in case 1 with faster finger flow is much smaller than that in case 2 with a slower volume-averaged manner. The results are consistent with the conclusion of Bianchi and Pedretti [43] that fingering is more pronounced when the relative entropy HR is low. Therefore, in a way similar to the derivation of geological entropy indicators, it can be recognized that the relative entropy HR, which takes into account the spatial distributions of the heterogeneous fracture network including connectivity, extent and spacing, and aperture distribution, has significant effects on the development and propagation of preferential flow during heterogeneous seepage analysis. Note that the relative entropy HR cannot be simply treated as an alternative to connectivity, but rather a comprehensive description including the connectivity, the spatial distribution of aperture, and extent and orientation.Figure 10
Variation ofHR with different square grid size.


## 5.4. Influences of Fracture Orientation on Rainfall Infiltration

The fourth example is to investigate the effect of rainfall infiltration on a fractured rock slope with two sets of random distribution fractures, where the rapid generation of higher hydraulic pressure within a short possible time may become one of the critical factors on slope stability. Four types of discrete fracture networks are employed and their geometrical parameters are given in Table3. The corresponding computational meshes are created using the Monte-Carlo method as presented in Figure 11. Suppose that the strong rainfall process lasts for 50 hours with a constant intensity of 0.001 m/hr. In the analysis, the time increment of each time step is adopted as 1.0 hr.Table 3
Parameters of fractures and probability model.
Dip (°)
Aperture (10−4 m)
Mean trace length (m)
Mean spacing (m)
Mean values
Variance
Probability model
Case1
0.1
Normal
0.2
15
0.3
2.0
135
0.1
0.8
Case2
30
0.3
2.0
135
0.1
0.8
Case3
45
0.3
2.0
135
0.1
0.8
Case4
60
0.3
2.0
135
0.1
0.8Figure 11
Numerical realizations of two fracture systems in a fractured rock slope under rainfall condition: (a) dip of 15 degrees; (b) dip of 30 degrees; (c) dip of 45 degrees; (d) dip of 60 degrees.
(a)
Case1
(b)
Case2
(c)
Case3
(d)
Case4The left and bottom boundaries are impermeable and the right boundary under the initial water table is always treated as a water pressure head boundary. The infiltration boundaries are imposed on the other boundaries along the slope. Due to the high gradient of the slope, it is reasonable to assume that no runoff would be produced on the infiltration boundary, which creates a zero value ofψpond. As an initial condition, we assume that the domain is under a hydrostatic state with an external water table situated at 5.00 m. The water pressure head below the groundwater level increases linearly from 0 m at the water table position to 5 m at the bottom, and that above the groundwater level decreases linearly from 0 m at the water table position to -20 m at the top of the slope.The normalized flow rate distributions inside the rock slope after 10 and 20 hours of water infiltration are shown in Figure12. The seepage flow occurs prevalently along bedding plane fractures caused by higher hydraulic conductivity of larger facture aperture; a portion of factures are still not active during rainfall infiltration. It can be observed that, with the increase of fracture dip, the propagation of water seepage extends inward more rapidly and the expansion of transient saturated region is driven faster into a deeper domain. The fractures with a steep dip can give rise to fast movement of water infiltration. A comparison of the normalized flow rate distributions ranging from 10 to 20 hr indicates that the orientation of fractures is also an important feature for describing the evolution of fast seepage flow. Figure 13 shows the variation of the total net flux into the unsaturated domain with time. When the fracture dip increases from 15 to 45 degrees, the growth of the total net flux increases substantially; conversely, the growth of the total net flux turns out to be smaller when the fracture dip reaches 60 degrees, especially after 20 hours water infiltration.Figure 12
The normalized flow rate distributions in the slope at times (a) 10 hr and (b) 20 hr.
(a)
Case1
(b)
Case2
(c)
Case3
(d)
Case4Figure 13
Evolution of the total net flow into the unsaturated domain.


## 6. Conclusion

As the permeability of matrix is extremely low and can be neglected, the fractures would have a great impact on the partially saturated water seepage flow. In order to understand the water seepage behavior in such heterogeneous fracture systems, a numerical approach based on the PVI formulations of variably saturated water flow in a discrete fracture network has been proposed.In consideration of the essential distinction in geometrical characteristic between the two-dimensional fracture and three-dimensional porous media, the constitutive relationships including water pressure head and relative permeability as functions of saturation, compared with the conventional van Genuchten’s model, are employed. While modeling of unsaturated flow in complex fracture systems is difficult and uncertain, the seepage surface and infiltration boundaries are unified as a complementary condition of Signorini-type formulation. Through the equivalence proof between the PDE and PVI formulations of the variably saturated water flow seepage problems in the fracture networks, the difficulty in solving such problems with boundary nonlinearity is reduced. In addition, the corresponding numerical finite element method is presented in detail.During the analysis of saturated-unsaturated water seepage simulation, the validity of the proposed procedure is demonstrated by the match between the predicted results of the proposed PVI method and analytical/numerical solutions. The model calculations in complicated fracture systems also suggest that water seepage may proceed by means of fast preferential flow paths along partial fractures due to the inhomogeneity of spatial distribution, which has a great disparity with the volume-averaged mode. Simultaneously, the variably saturated seepage behavior is significantly sensitive to the geometrical characteristics of the fracture system. The fast flow pathways are not only dominated by the spatial distribution of fracture aperture, but also strongly dependent on the connectivity, density, and orientation. Actually, the water flow through complicated fracture networks is also affected by the fracture density and connectivity. In particular, the relative entropyHR is adopted to describe the heterogeneity in the fracture networks, and it is found that fingering is more pronounced when the relative entropy HR is lower. However, more investigations are still required to evaluate the dependence of the relative entropy HR on the saturated-unsaturated flow in fracture networks. Further efforts are needed to quantify the relationships between water flow behavior and fracture spatial geometry (including the distribution of aperture, extent, and orientation) in random or correlated (fractal) fracture networks. These aspects will be investigated in our future study.



---

*Source: 9062569-2018-01-11.xml*