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.
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 Section 2, 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 Figure 1, 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 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 , including the net gain from fluid sources. With this extension, the governing equation for an arbitrary fracture segment , as shown in Figure 2, can be expressed in a modified form of Richards’ equation: where and are the pressure head [L] and the elevation heads [L], respectively, in the fracture segment ; is the relative permeability; is the specific moisture capacity of fractures [L−1]; is the water saturation. The saturated hydraulic conductivity [LT−1] of an individual fracture segment follows the cubic law, which can be stated as follows:where denotes the gravity acceleration [LT−2] and is the dynamic viscosity of water [L2T]; is the hydraulic fracture aperture [L], rather than the mean geometric aperture of the fracture.
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 [L] is commonly adopted to replace the mean aperture in (2).
As shown in Figure 3, 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 fracture segments meeting at intersection i, the sum of the total inflow flux and outflow flux at intersection should be zero based on the principle of mass conservation, which can be expressed aswhere is an extension of Darcy flow velocity [LT−1] within the fracture segment under variable saturated flow condition, given by .
2.2. Constitutive Relationships
Constitutive relationships between the pressure head , relative permeability , and water saturation 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 where , n, and are empirical fitting parameters with , is the effective saturation given by , and and 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 , the relevant saturation and relative permeability yieldwhere the local aperture larger than cannot be occupied by the water flow, is the interfacial tension or surface tension, is the contact angle, 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 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 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 for unsaturated seepage flow analysis in a fracture network can be described according to
On the rain infiltration surfaces , 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 aswhere is the rainfall intensity [LT−1] and 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:where and refer to critical flow rate [L2T−1] and water pressure head [L], respectively. From the above unified expression, we have on the seepage boundary and on the infiltration boundary.
In addition, the saturated-unsaturated water flow in fracture networks should satisfy the following conditions:(1)The initial condition:(2)The pressure head boundary condition(3)The flux boundary condition
where is the prescribed water head [L] and is the prescribed flow rate [L2T−1] for fracture .
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 functionsthe theoretical solution with respect to should satisfywhere
By employing integration by parts, (15) is given bywhere is the total number of fracture nodes.
Inserting (13), (16), and (17) into (11), the following relationship can be obtained:
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 gives
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 . Inserting and into (18), where is an arbitrary function whose values are zeros at all the nodes in the whole fracture network, the continuity equation (1) can be derived:
Thus, (18) can be expressed as
Supposing and in (19), where is any function that becomes zero at the nodes on , , and , the mass conservation equation on the intersections can be obtained:
Equation (21) can be rewritten as
Taking and in (21), with being any function that equals zero at the nodes of boundaries and , the flux condition yields
Herein, (23) can be given as
As the term in the above inequality is a definite value, the inequality condition of (23) for leads to
From (12), is an arbitrary function that satisfies for , so the following equation is obtained:
Particularly, inserting in (23) results in
From (27) and (28) and the fact that for , the complementary condition of Signorini-type formulation for is expressed as
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 , 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 corresponding to in a trial vector spacewhich should obey the following relationships:wherewhere is the time step, is the time increment, and and 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: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 m−1, , and . The other two examples are based on the work of Ye et al. [16, 31] and the relevant parameters are and . The three related retention curves are presented in Figure 4.
(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 m pressure head.
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 Figure 5. 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
Example 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 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].
(a) Homogenous dam
(b) Nonhomogenous dam
Figure 6 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 Figure 7. 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.
(a) Case
(b) Case
Different 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 Table 1. 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 is ten times that in case , and the mean spacing in case is half that in case . The corresponding realizations of two different fracture systems are presented in Figure 7. For numerical modeling, a time step of 1 s is adopted.
Figure 8 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.
(a) Case
(b) Case
The 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 cases and , 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 yieldwhere is the relative entropy index, which is a quantitative measure of the spatial disorder in a discrete random system; is the global entropy of the discrete random system representing a set of grid blocks; is the local entropy of each square grid; ’s are defined as the volumetric proportions of the fractures or rock masses over the whole domain of interest; ’s are defined as the volumetric proportions of the fractures or rock masses over each square grid; is the average of the local entropy over the entire domain.
An example to characterize such fracture network system is shown in Figure 9; 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, ’s of fractures and rock masses are equal to 0.02 and 0.98 according to(35), corresponding to an value of 0.098. When the two systems are divided into nine subsets as shown in Figure 9, the relative entropy values are reported in Table 2. The index for the connected system is equal to 0.867. On the other hand, the unconnected system has a very large value, which is 0.969. It is found that the indicator is sensitive to the connectivity of the fracture networks. Furthermore, Bianchi and Pedretti [43] pointed out that is a global metric that integrates multiple properties of the spatial distribution of hydraulic conductivity.
(a) Connected ()
(b) Unconnected ()
Herein, the relative entropy 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, increases with the size of square grid. However, the value in case with faster finger flow is much smaller than that in case 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 is low. Therefore, in a way similar to the derivation of geological entropy indicators, it can be recognized that the relative entropy , 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 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.
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 Table 3. 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.
(a) Case
(b) Case
(c) Case
(d) Case
The 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 . 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 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 Figure 12. 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.
(a) Case
(b) Case
(c) Case
(d) Case
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 entropy is adopted to describe the heterogeneity in the fracture networks, and it is found that fingering is more pronounced when the relative entropy is lower. However, more investigations are still required to evaluate the dependence of the relative entropy 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.
Nomenclature
| Pressure head, L | |
| Elevation heads, L | |
| Gravity acceleration, LT−2 | |
| Darcy flow velocity, LT−1 | |
| Flux rate, L2T−1 | |
| Effective saturation | |
| Residual saturation | |
| Satiated saturation | |
| Shape function | |
| Iteration step | |
| Error tolerance | |
| Time step | |
| Time increment, T | |
| Hydraulic aperture of fracture segment , L | |
| Relative permeability | |
| Water saturation | |
| Rainfall intensity, LT−1 | |
| Ponding depth, L | |
| Critical flow rate, L2T−1 | |
| Prescribed water head, L | |
| Prescribed flow rate, L2T−1 | |
| Saturated hydraulic conductivity, LT−1 | |
| Dynamic viscosity of water, L2T | |
| Equivalent hydraulic aperture, L | |
| Critical water pressure head, L | |
| Specific moisture capacity, L−1. |
Disclosure
All data used for the plots were generated using the authors’ own computer codes. The data are available upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
The financial supports from the National Natural Science Foundation of China (nos. 51709207, 41762020, 51679173, and 51604195) and Natural Science Foundation of Hubei Province (no. 2015CFA142) are gratefully acknowledged.