# A Novel Numerical Model for Fluid Flow in 3D Fractured Porous Media Based on an Equivalent Matrix-Fracture Network

Authors: Chi Yao; Chen He; Jianhua Yang; Qinghui Jiang; Jinsong Huang; Chuangbing Zhou
Journal: Geofluids
 (2019)

---

## Abstract

An original 3D numerical approach for fluid flow in fractured porous media is proposed. The whole research domain is discretized by the Delaunay tetrahedron based on the concept of node saturation. Tetrahedral blocks are impermeable, and fluid only flows through the interconnected interfaces between blocks. Fractures and the porous matrix are replaced by the triangular interface network, which is the so-called equivalent matrix-fracture network (EMFN). In this way, the three-dimensional seepage problem becomes a two-dimensional problem. The finite element method is used to solve the steady-state flow problem. The big finding is that the ratio of the macroconductivity of the whole interface network to the local conductivity of an interface is linearly related to the cubic root of the number of nodes used for mesh generation. A formula is presented to describe this relationship. With this formula, we can make sure that the EMFN produces the same macroscopic hydraulic conductivity as the intact rock. The approach is applied in a series of numerical tests to demonstrate its efficiency. Effects of the hydraulic aperture of fracture and connectivity of the fracture network on the effective hydraulic conductivity of fractured rock masses are systematically investigated.


---

## Body

## 1. Introduction

The research for seepage flow in the subsurface system is of great significance in many geological and environmental engineering applications, for example, petroleum and gas exploitation, greenhouse gas storage, and nuclear waste treatment. Most rocks in nature are cut by fracture or fracture network, and the fractures have significant influence on the deformation and seepage characteristics of rocks [1–4]. It is well known that fluid tends to flow through the easiest path to the downstream. For fractured rock, the conductivity of fractures is often several orders of magnitude higher than that of the porous matrix, which means that fractures and fracture networks can strongly influence or even dominate the conductivity and flow paths of rock masses. Moreover, fluid exchange between the porous matrix and fractures often results in complex hydraulic behaviors, such as seepage flow in aquifers and reservoirs. It is generally hard to solve seepage problems in such complicated materials with analytical methods. On the other hand, many numerical simulation methods, such as finite element method [5, 6], finite volume method [7, 8], boundary element method [9], and numerical manifold method [10, 11], provide efficient tools for the analysis of fluid flow in fractured porous media and can help us better predict the hydraulic properties of underground engineering.There are three major kinds of model for fluid flow simulation in fractured porous media which based on equivalent continuum/porous media (EPM), discrete fracture network (DFN), or fractured porous media (FPM). The EPM-based model assumes that fluid flow obeys Darcy’s law, in which fractures and matrix are approximated as a continuum and the classical continuous seepage theory is used for analysis [12–14]. This kind of model provides a phenomenological description of flow in fractured porous media in large dimension. However, it is limited by whether the REV exists and cannot reflect exactly the local conductivity anisotropy induced by fractures and the preferential flow effect within fracture networks. DFN-based models assume that rock blocks are impermeable and flow movement only takes place in fracture networks [15–17]. The geometry and topology of fracture networks are generally described explicitly by this kind of approach. DFN-based models are applicable only when the fracture network is well-connected, and the conductivity of the rock matrix is negligible. Methods based on the FPM concept consider that fluid flows simultaneously within fractures and the rock matrix [7, 11, 18–21]. This kind of approach inherits the advantages of EPM and DFN. It not only describes the geometry of the fracture network but also takes into consideration the contribution from the matrix and the fluid exchange between matrix and fractures. However, this method has difficulties in numerical mesh generation due to geometry complexity of fracture networks, and the computational cost is generally very high. The three kinds of model mentioned previously have their advantages in fluid seepage simulation, respectively, among which the FPM-based models are usually more comprehensive, more consistent with physical reality, and thus could provide more details.To overcome the drawbacks of three previously mentioned approaches, Yao et al. [22] proposed a novel equivalent matrix-fracture network (EMFN) approach, in which flows in the rock matrix and fractures are represented by the flow in an equivalent fracture network. However, Yao et al. [22] only considered 2D problems where the Voronoi diagram was used for discrete 2D domains and the one-dimensional equivalent fracture network was used to study the conductivity evolution of fractured rock masses. It is well known that the 2D fracture network tends to underestimate the conductivity compared with those based on 3D description since the connectivity of the 2D model is inferior to the 3D model [5, 23, 24]. As a result, 3D models are necessary for the prediction of hydraulic properties of fractured porous media.In this study, the EMFN approach is further developed to model 3D fluid flows in fractured porous media. The porous matrix and fracture network are replaced by an equivalent matrix-fracture network (EMFN). Delaunay tetrahedrons are used to construct the 3D fracture network. It is shown that Delaunay tetrahedrons are efficient and able to handle irregular boundary geometry and complex fracture distributions. The proposed approach is applied to study the conductivity of complex fractured rock masses. Results show that there is a unique relationship between the macroscopic conductivity of the whole interface network and the microscopic conductivity of interfaces, the ratio of which is linearly related to the cubic root of the number of nodes used for mesh generation. With this relationship, the so-called equivalent matrix-fracture network (EMFN) can produce the same conductivity of intact rock, which lays the foundation for the proposed model. The effectiveness of this approach is then demonstrated by several examples.


## 2. The Equivalent Matrix-Fracture Network Model

### 2.1. The Numerical Model

In this approach, rock masses and fractures are discretized into a triangular interface network in the 3D domain. Fluid flow occurs on the network and follows Darcy’s law. Firstly, a 2D model is used to illustrate the concept of the equivalent matrix-fracture network model. A 2D schematic of one streamline (solid red line) in the model is depicted in Figure1(a), where the triangular interfaces are reduced to line segments. A 3D conceptual diagram of the interface network is also shown in Figure 1(b). The triangular interfaces between tetrahedrons are the paths through which fluid flows. Rock matrix and fractures are represented by interface networks in the same mesh. The fracture thickness is not considered. The only difference between the matrix interface and the fracture interface is that they are assigned different transmissivities while building the equivalent equation. In this way, the pore structure of the porous media and the fracture network is replaced by an equivalent interface network, which is called the 3D equivalent matrix-fracture network (EMFN). The core concept of the EMFN model is to use a sufficiently complex and chaotic discrete and permeable interface network to fill a three-dimensional volume, making the effective permeability of the interface network equivalent to the continuous media. Then, the fractures and matrix are simultaneously represented as lower-dimensional domains (2D surfaces) in the volume, and the Galerkin-based finite element method is used to solve the flow problem.Figure 1
Concept illustration of fluid flow in an equivalent matrix-fracture network model in 2D and 3D.
(a)
(b)Assuming that the transmissivity of an individual interface isTe, the continuity equation of the triangular interface abc under steady-state flow conditions in the local coordinate (Figure 2(a)) is expressed as
(1)Te∂2h∂x′2+∂2h∂y′2=0,where h is the hydraulic head, 1 m of which is equivalent to 9.81 kPa. Suppose an intersection line is the common edge of two triangular fracture surfaces (Figure 2(b)). Neglecting the storage capacity of the porous media, for any node on this line, the equation of mass conservation is
(2)∑i=1mqni=0,where m is the number of interfaces passing through this node, and qni is the flow rate from this node to the interface i,
(3)qni=−T∂h∂ni.Figure 2
The local illustration of triangular elements.
(a)
(b)Illustrated in Figure2(b), ni is the direction vector that is perpendicular to the intersection line and parallel to the plane i. The flow problem in the whole domain also needs to satisfy the hydraulic head boundary conditions and the flux boundary conditions.The finite element method is used to solve the steady-state flow problem in the EMFN model. The hydraulic headh of any node on a triangular element can be represented by the hydraulic head at three vertices (named i,j,k) surrounding this element:
(4)h=Nihi+Njhj+Nkhk,where Ni, Nj, and Nk are the shape functions. After the finite element approximation, the total equivalence equation of the steady-state flow is built up:
(5)KH=F,where
(6)K=∑ΩKe,Ke=∬eBTTeBdx′dy′,B=∂Ni∂x′∂Nj∂x′∂Nk∂x′∂Ni∂y′∂Nj∂y′∂Nk∂y′,where x′, y′, and z′ are local coordinates and H and F are, respectively, the unknown hydraulic head vector and the vector calculated from the known boundary condition. The solution of the whole seepage flow field can be obtained by solving (5).


### 2.2. Three-Dimensional Mesh Generator for Complex Fracture Network Based on Delaunay Tessellation

In the meshing process of the research domain with irregular boundary geometry and incorporating complex fracture distribution, auto-unstructured mesh generation techniques, such as Delaunay triangulation (in 2D) and tetrahedralization (in 3D) [7, 25], are relatively appropriate and efficient. In this paper, the 3D research domain is discretized by Delaunay tetrahedrons (Figure 1(b)). Delaunay tetrahedralization is a fast and relatively high-quality procedure for automatic spatial subdivision. Once the spatial node set and the geometric boundary condition have been given, the unique tetrahedral mesh can be generated. The node set exists as vertexes of tetrahedrons in the mesh.A random node saturation procedure is introduced in the node set generation illustrated in Figure3. A parameter lmin is used to define the minimum distance between any two nodes in the mesh, then the node keeps being put one by one randomly into the domain until no other node meets the condition. Then, the open-source Delaunay tetrahedralization code, TetGen [26], is applied to discretize the domain based on the node set. The fracture geometry can be described by a special Delaunay constrained facet and edge generation strategy. The process of the mesh generation program is as follows:
Step 1:
Define a research domain (in this paper, the domain is always a cube with side lengthL=100m as shown in Figure 1(b)), and then nodes with a random 3D coordinate were put into this domain one by one. Some of the nodes are put in firstly because they are used to form the fracture planes, intersection lines, and domain boundaries; the others are used to form the rock matrix and filled the rest of the domain.
Step 2:
Generate a strict Delaunay tetrahedral mesh using the existing node set.
Step 3:
Check the integrity of the fracture planes. If it is found that the area of some planes is missing, call the local node infill algorithm (described later in this section) and jump back to the Step 2.
Step 4:
Output the mesh data.Figure 3
Illustration of the random node saturation procedure.It needs to be pointed out that, in Step 1, almost all node insertion processes follow the node saturation principle. However, due to the presence of intersecting cracks with a very small dihedral angle, the nodes near their intersection line may break this principle and have a smallerlmin. Fortunately, the number of abnormal nodes is far smaller than the total number of nodes when the fracture network is not quite dense, which means that the effect of the abnormal nodes on the calculation of rock conductivity is negligibly small. This problem needs to be avoided through a graded-mesh-based approach in the future.Here, we introduce the local node infill algorithm in Step 3. For a node set (including the fracture geometry), the mesh generated by Delaunay’s method is unique (the situation of 5 or more nodes on one spherical surface can be avoided by disturbing the coordinates of nodes). Unfortunately, some of the tetrahedrons adopted may pass through these fracture planes, resulting in the loss of the area of some fracture planes. A schematic fracture plane is illustrated in Figure4(a), which lost a triangular area. The reason is shown in Figure 4(b): three tetrahedrons (pqab, pqbc, and pqca) passed through the plane, causing the plane to lose a triangular area abc, where pq is the common edge of three tetrahedrons.Figure 4
Illustration of a local node infill process.
(a)
(b)
(c)In rigid Delaunay tetrahedralization, the circumscribed sphere of any tetrahedron includes no other nodes besides four vertexes of the tetrahedron. So, one can appropriately put some nodes into the domain appropriately to break these undesirable tetrahedrons and re-form new tetrahedrons that do not pass through fracture planes. In fact, these nodes are the intersection nodes of fracture planes with the common edge of abnormal tetrahedrons, for example, the intersection node ofpq with the plane abc in Figure 4(c). The final mesh is illustrated in Figure 5, where the green interface represents the fracture planes and the blue interface represents the equivalent matrix (for clarity, some interfaces close to the boundary are neglected).Figure 5
The final mesh.


### 2.3. Relationship between Mesh Parameter and Macroscopic Conductivity of EDFN Model

By adopting the random node saturation method, the density of the node set is determined throughlmin. In fact, there is a linear relationship between lmin and the total node number N of the cubic domain with side length L:
(7)N3≈0.8649Llmin+1.0180.Formula (7) is derived from 10 groups of numerical tests, in which L/lmin=10~200, and each group contains 8 specimens. Numerical results are shown in Figure 6; the solid line is the result of a linear fit. Moreover, the coefficient of variation of N of one group is, respectively, 1.10%, 0.12%, and 0.01%, when L/lmin equals 10, 50, and 200.Figure 6
Relationship betweenL/lmin and N3.In the EMFN model, the pore system of the rock matrix is replaced by an interconnected equivalent matrix-fracture network, which is constructed by an assembly of triangular interfaces. The hydraulic conductivityTm of each equivalent matrix-fracture is used to ensure this system is hydraulically equivalent to the rock matrix. In this paper, all Tm are assumed to be identical. To verify the effectiveness of this method and to study the mesh dependency of the EMFN model at the same time, a group of numerical experiments is designed. There are 6 groups (each group contains 8 specimens) of meshes without any fracture, and the only difference between groups is lmin. The research domain is cubic and L=100m. A pair of opposite boundaries are hydraulic head boundaries, and the others are impervious boundaries. According to Darcy’s law,
(8)Q=AkJ,Here,J=h1−h2/L=1m/m. Defining dimensionless quantity θ to evaluate the transmission capacity of the research domain relating to Tm,
(9)θi=QiJLTm,where the subscript i represents the i direction (x-, y-, or z-direction) in the global coordinate system. The results are listed in Table 1, where N is the average number of nodes, θi is the average value of one group in the i direction, and C.V.θi is the coefficient of variation of θx, θy, and θz in a group. One can see that the Delaunay mesh generated from a randomly and uniformly distributed set of nodes with identical Tm on each equivalent fracture has almost the same θ in three dimensions, which means the equivalent permeability of the matrix-fracture network generated by the process in Section 2.2 is almost isotropic. Further, the variation coefficient of conductivity becomes smaller with the decrease in lmin.Table 1
Statistical analysis of conductivity for different numerical tests.
#
lmin
N
θx
θy
θz
C.V.θi
1
22.18
118
21.79
21.95
21.74
0.498%
2
9.55
1011
47.06
47.13
47.19
0.133%
3
3.01
26,426
145.24
145.31
145.21
0.036%
4
1.78
122,453
244.10
244.15
244.08
0.016%
5
1.26
335,968
342.83
342.81
342.75
0.012%
6
0.98
716,968
442.03
442.04
442.04
0.002%From Figure7, one can see that θ and N3 have a linear relationship:
(10)θ=4.9273N3.Figure 7
Relationship betweenθ and N3.Substituting (10) into (9) results in
(11)Q=4.9273JLTmN3.According to Darcy’s theorem, macroscopic flow rateQ and macroscopic conductivity coefficient k have the following relationship:
(12)Q=kJL2.It can be derived from (11) and (12) that
(13)k=4.9273TmNL33.When using the mesh density parameterlmin as an argument, (13) can be rewritten as
(14)k=4.9273TmL0.8649Llmin+1.0180.At present, the relationship between the macroscopic conductivity coefficient and the microscopic parameter of the EMFN model has been established.


### 2.4. A Modified Method for Complex Fracture Network

The functionkTm,L,lmin describes the relationship between macroscopic hydraulic parameters and the microscopic parameters in EMFN. It is obtained from data fitting for the mesh that does not include any fracture, which means that the nodes used to generate the Delaunay-based mesh satisfy the “random node saturation procedure” in any part of the research domain. In other words, the density of chaos-oriented interfaces is homogeneous, and the influence of a slight disturbance in node density is negligible. However, when it comes to complex fracture networks (such as those shown in Figure 8, where the green disks are fractures and the blue cube is the research domain), the effect of additional refinement nodes cannot be ignored because they are too many. Where the refinement nodes are added, the density of the interface will be higher than that of other areas where they do not contain refinement nodes. Refinement nodes exist only in the fracture plane. Under this feature, according to the formation rule of Delaunay’s tetrahedron and the EMFN model, the equivalent permeability of the matrix will be affected only within the distance of one tetrahedron around the refinement node, and the transmissivity within the fracture plane will not be affected. The denser the interface network is, the higher the matrix network equivalent permeability is.Figure 8
Illustration for a complex fracture network (Nfr=10, 40, 80).
(a)
(b)
(c)In light of the problem, a practical strategy to reduce the influence of refinement nodes on the equivalent permeability of the matrix network is developed. We mark all the refinement nodes added by Step 3 introduced in Section2.2 and then find out in all the matrix interfaces that at least one of their vertex nodes is the refinement node. Finally, the conductivity of these interfaces is reduced because they are denser than we expected. Through a series of tests, a pragmatic reduction coefficient, 0.55, was determined. That is, if Tm is the transmissivity of matrix interfaces used in (14), then the transmissivity of these dense matrix interfaces is 0.45 Tm. We cannot prove mathematically that the reduction of transmissivity is correct, but it works.A series of experiments were carried out to verify the effectiveness of the reduction strategy. We designed a research domain thatlmin/L=0.02, and in which the centers of fractures are randomly distributed, and their orientation is random and isotropic with periodic boundary conditions (i.e., the example in Figure 8). The radius of fractures is identical and equal to 0.25 L. Under these conditions, samples with different fracture densities were tested. Fracture densities are assessed using the number Nfr of fractures inside the research domain. In this test, all the transmissivities of interfaces are identical and equal to Tm. Thus, although fractures are indeed described in the mesh, the fracture interfaces are no different from the matrix interfaces. So the influence of nonhomogeneous interface density and the effect of transmissivity reduction strategy can be easily known through steady-state flow simulation. We define the hydraulic conductivity of the matrix calculated by (14) as K0 and the hydraulic conductivity of fractured media obtained through EMFN as K. As shown in Figure 9, the black squares are K/K0 under the reduction strategy, which does not change significantly with the increase in fracture number (more specifically, the rise of refinement node). However, the results (red spots) obtained from the samples that not have a reduction procedure are obviously larger. The percentages of refinement nodes in total nodes are also noted in Figure 9.Figure 9
Effects of reduction strategy for samples with differentNfr.


### 2.5. Verification of the EMFN Model

In the EMFN model, either fracture interfaces or the equivalent matrix interfaces are all treated in the same way of the equivalent-fracture network, in which the “fracture interfaces” and the “matrix interfaces” are givenT independently; for fracture it is Tf, and for matrix it is Tm. Unfortunately, for flow problems in the discrete fracture porous media, it is hard to find an analytic solution of conductivity even for very simple conditions [27]. However, we can use the superposition principle in conductivity estimation to estimate the conductivity in some well-designed conditions. Two examples are given in this section. The first one is to test the conductivity of the rock a single rectangular fracture passing through the domain; the other one is to examine the conductivity tensor of the rock which has a disk fracture embedded with varying orientations in the centroid of the domain.When a rectangular fracture is parallel to the hydraulic gradientJ (sketching is depicted in Figure 10), its contribution to flow rate in the direction of J is easy to get by
(15)Qf=TflfJ,where Qf is the flow rate through fracture and lf is the trace length of the fracture. In this example, lf=60m, Tf=1×10−8m2/s, Tm=1×10−12m2/s, and L/lmin=50. Substituting (15) into (9) results in the x-direction as θfx=6×103. The transmission capacity of the matrix can be calculated by (9), (12), and (14) as θmx=2.18×102. The result of simulation using the EMFN model is θx=6.217×103.Figure 10
Single rectangular crack throughout the domain.
(a)
(b)The second example is illustrated in Figure11(a): a disk fracture (30 m in diameter) is embedded in the cubic domain, Tf=1×10−8m2/s, Tm=1×10−12m2/s, and L/lmin=50. The boundary condition is depicted in Figure 11(c), and a pair of boundaries perpendicular to the x-axis are given, respectively, constant hydraulic head h1=100m and h2=0m; others are corresponding linear variable hydraulic head boundaries. The center of the disk fracture coincides with the centroid of the cubic domain; the normal vector n of the fracture is always perpendicular to the z-axis. α is defined as the angle formed by the x-axis and n. The conductivity tensors when α equal to 0°, 15°, 30°, 45°, 60°, 75°, and 90° are calculated.Figure 11
Single circular fracture rotating in the center.
(a)
(b)
(c)The conductivity tensorK is the function of α and can be expressed as
(16)Kα=KxxαKxyα0KyxαKyyα000Kzzα.Whenα=0°,
(17)Kα=K1000K2000K3,where K1,K2,andK3 are, respectively, the coefficient of conductivity in the x,y,andz directions. Components in (16) can be written as
(18)Kxxα=K1cos2α+K2sin2α,Kyyα=K1sin2α+K2cos2α,Kxyα=Kyxα=K1−K2sinαcosα,Kzzα=K3.SubstitutingK1,K2,andK3 adopted from the numerical simulation when α=0° (here K1=2.111×10−12m/s, K2=K3=2.144×10−12m/s) into (16), (17), and (18), the relationships between α and components of K are received. K0 is defined as the coefficient of conductivity simulated in the same conditions but without any fracture (in this case, K0=2.111×10−12m/s and equal to K1, which means the fracture has no improvement on the conductivity in the direction parallel to n). The comparison between the numerical result using the EMFN model and the analytical result calculated by (16), (17), and (18) is plotted in Figure 12. One can see that the numerical results obtained are very close to the theoretical prediction values, which proves the validity of the EMFN method proposed in this paper.Figure 12
The comparison theoretical (solid line) and the numerical (symbol) result of components ofK.


## 2.1. The Numerical Model

In this approach, rock masses and fractures are discretized into a triangular interface network in the 3D domain. Fluid flow occurs on the network and follows Darcy’s law. Firstly, a 2D model is used to illustrate the concept of the equivalent matrix-fracture network model. A 2D schematic of one streamline (solid red line) in the model is depicted in Figure1(a), where the triangular interfaces are reduced to line segments. A 3D conceptual diagram of the interface network is also shown in Figure 1(b). The triangular interfaces between tetrahedrons are the paths through which fluid flows. Rock matrix and fractures are represented by interface networks in the same mesh. The fracture thickness is not considered. The only difference between the matrix interface and the fracture interface is that they are assigned different transmissivities while building the equivalent equation. In this way, the pore structure of the porous media and the fracture network is replaced by an equivalent interface network, which is called the 3D equivalent matrix-fracture network (EMFN). The core concept of the EMFN model is to use a sufficiently complex and chaotic discrete and permeable interface network to fill a three-dimensional volume, making the effective permeability of the interface network equivalent to the continuous media. Then, the fractures and matrix are simultaneously represented as lower-dimensional domains (2D surfaces) in the volume, and the Galerkin-based finite element method is used to solve the flow problem.Figure 1
Concept illustration of fluid flow in an equivalent matrix-fracture network model in 2D and 3D.
(a)
(b)Assuming that the transmissivity of an individual interface isTe, the continuity equation of the triangular interface abc under steady-state flow conditions in the local coordinate (Figure 2(a)) is expressed as
(1)Te∂2h∂x′2+∂2h∂y′2=0,where h is the hydraulic head, 1 m of which is equivalent to 9.81 kPa. Suppose an intersection line is the common edge of two triangular fracture surfaces (Figure 2(b)). Neglecting the storage capacity of the porous media, for any node on this line, the equation of mass conservation is
(2)∑i=1mqni=0,where m is the number of interfaces passing through this node, and qni is the flow rate from this node to the interface i,
(3)qni=−T∂h∂ni.Figure 2
The local illustration of triangular elements.
(a)
(b)Illustrated in Figure2(b), ni is the direction vector that is perpendicular to the intersection line and parallel to the plane i. The flow problem in the whole domain also needs to satisfy the hydraulic head boundary conditions and the flux boundary conditions.The finite element method is used to solve the steady-state flow problem in the EMFN model. The hydraulic headh of any node on a triangular element can be represented by the hydraulic head at three vertices (named i,j,k) surrounding this element:
(4)h=Nihi+Njhj+Nkhk,where Ni, Nj, and Nk are the shape functions. After the finite element approximation, the total equivalence equation of the steady-state flow is built up:
(5)KH=F,where
(6)K=∑ΩKe,Ke=∬eBTTeBdx′dy′,B=∂Ni∂x′∂Nj∂x′∂Nk∂x′∂Ni∂y′∂Nj∂y′∂Nk∂y′,where x′, y′, and z′ are local coordinates and H and F are, respectively, the unknown hydraulic head vector and the vector calculated from the known boundary condition. The solution of the whole seepage flow field can be obtained by solving (5).


## 2.2. Three-Dimensional Mesh Generator for Complex Fracture Network Based on Delaunay Tessellation

In the meshing process of the research domain with irregular boundary geometry and incorporating complex fracture distribution, auto-unstructured mesh generation techniques, such as Delaunay triangulation (in 2D) and tetrahedralization (in 3D) [7, 25], are relatively appropriate and efficient. In this paper, the 3D research domain is discretized by Delaunay tetrahedrons (Figure 1(b)). Delaunay tetrahedralization is a fast and relatively high-quality procedure for automatic spatial subdivision. Once the spatial node set and the geometric boundary condition have been given, the unique tetrahedral mesh can be generated. The node set exists as vertexes of tetrahedrons in the mesh.A random node saturation procedure is introduced in the node set generation illustrated in Figure3. A parameter lmin is used to define the minimum distance between any two nodes in the mesh, then the node keeps being put one by one randomly into the domain until no other node meets the condition. Then, the open-source Delaunay tetrahedralization code, TetGen [26], is applied to discretize the domain based on the node set. The fracture geometry can be described by a special Delaunay constrained facet and edge generation strategy. The process of the mesh generation program is as follows:
Step 1:
Define a research domain (in this paper, the domain is always a cube with side lengthL=100m as shown in Figure 1(b)), and then nodes with a random 3D coordinate were put into this domain one by one. Some of the nodes are put in firstly because they are used to form the fracture planes, intersection lines, and domain boundaries; the others are used to form the rock matrix and filled the rest of the domain.
Step 2:
Generate a strict Delaunay tetrahedral mesh using the existing node set.
Step 3:
Check the integrity of the fracture planes. If it is found that the area of some planes is missing, call the local node infill algorithm (described later in this section) and jump back to the Step 2.
Step 4:
Output the mesh data.Figure 3
Illustration of the random node saturation procedure.It needs to be pointed out that, in Step 1, almost all node insertion processes follow the node saturation principle. However, due to the presence of intersecting cracks with a very small dihedral angle, the nodes near their intersection line may break this principle and have a smallerlmin. Fortunately, the number of abnormal nodes is far smaller than the total number of nodes when the fracture network is not quite dense, which means that the effect of the abnormal nodes on the calculation of rock conductivity is negligibly small. This problem needs to be avoided through a graded-mesh-based approach in the future.Here, we introduce the local node infill algorithm in Step 3. For a node set (including the fracture geometry), the mesh generated by Delaunay’s method is unique (the situation of 5 or more nodes on one spherical surface can be avoided by disturbing the coordinates of nodes). Unfortunately, some of the tetrahedrons adopted may pass through these fracture planes, resulting in the loss of the area of some fracture planes. A schematic fracture plane is illustrated in Figure4(a), which lost a triangular area. The reason is shown in Figure 4(b): three tetrahedrons (pqab, pqbc, and pqca) passed through the plane, causing the plane to lose a triangular area abc, where pq is the common edge of three tetrahedrons.Figure 4
Illustration of a local node infill process.
(a)
(b)
(c)In rigid Delaunay tetrahedralization, the circumscribed sphere of any tetrahedron includes no other nodes besides four vertexes of the tetrahedron. So, one can appropriately put some nodes into the domain appropriately to break these undesirable tetrahedrons and re-form new tetrahedrons that do not pass through fracture planes. In fact, these nodes are the intersection nodes of fracture planes with the common edge of abnormal tetrahedrons, for example, the intersection node ofpq with the plane abc in Figure 4(c). The final mesh is illustrated in Figure 5, where the green interface represents the fracture planes and the blue interface represents the equivalent matrix (for clarity, some interfaces close to the boundary are neglected).Figure 5
The final mesh.


## 2.3. Relationship between Mesh Parameter and Macroscopic Conductivity of EDFN Model

By adopting the random node saturation method, the density of the node set is determined throughlmin. In fact, there is a linear relationship between lmin and the total node number N of the cubic domain with side length L:
(7)N3≈0.8649Llmin+1.0180.Formula (7) is derived from 10 groups of numerical tests, in which L/lmin=10~200, and each group contains 8 specimens. Numerical results are shown in Figure 6; the solid line is the result of a linear fit. Moreover, the coefficient of variation of N of one group is, respectively, 1.10%, 0.12%, and 0.01%, when L/lmin equals 10, 50, and 200.Figure 6
Relationship betweenL/lmin and N3.In the EMFN model, the pore system of the rock matrix is replaced by an interconnected equivalent matrix-fracture network, which is constructed by an assembly of triangular interfaces. The hydraulic conductivityTm of each equivalent matrix-fracture is used to ensure this system is hydraulically equivalent to the rock matrix. In this paper, all Tm are assumed to be identical. To verify the effectiveness of this method and to study the mesh dependency of the EMFN model at the same time, a group of numerical experiments is designed. There are 6 groups (each group contains 8 specimens) of meshes without any fracture, and the only difference between groups is lmin. The research domain is cubic and L=100m. A pair of opposite boundaries are hydraulic head boundaries, and the others are impervious boundaries. According to Darcy’s law,
(8)Q=AkJ,Here,J=h1−h2/L=1m/m. Defining dimensionless quantity θ to evaluate the transmission capacity of the research domain relating to Tm,
(9)θi=QiJLTm,where the subscript i represents the i direction (x-, y-, or z-direction) in the global coordinate system. The results are listed in Table 1, where N is the average number of nodes, θi is the average value of one group in the i direction, and C.V.θi is the coefficient of variation of θx, θy, and θz in a group. One can see that the Delaunay mesh generated from a randomly and uniformly distributed set of nodes with identical Tm on each equivalent fracture has almost the same θ in three dimensions, which means the equivalent permeability of the matrix-fracture network generated by the process in Section 2.2 is almost isotropic. Further, the variation coefficient of conductivity becomes smaller with the decrease in lmin.Table 1
Statistical analysis of conductivity for different numerical tests.
#
lmin
N
θx
θy
θz
C.V.θi
1
22.18
118
21.79
21.95
21.74
0.498%
2
9.55
1011
47.06
47.13
47.19
0.133%
3
3.01
26,426
145.24
145.31
145.21
0.036%
4
1.78
122,453
244.10
244.15
244.08
0.016%
5
1.26
335,968
342.83
342.81
342.75
0.012%
6
0.98
716,968
442.03
442.04
442.04
0.002%From Figure7, one can see that θ and N3 have a linear relationship:
(10)θ=4.9273N3.Figure 7
Relationship betweenθ and N3.Substituting (10) into (9) results in
(11)Q=4.9273JLTmN3.According to Darcy’s theorem, macroscopic flow rateQ and macroscopic conductivity coefficient k have the following relationship:
(12)Q=kJL2.It can be derived from (11) and (12) that
(13)k=4.9273TmNL33.When using the mesh density parameterlmin as an argument, (13) can be rewritten as
(14)k=4.9273TmL0.8649Llmin+1.0180.At present, the relationship between the macroscopic conductivity coefficient and the microscopic parameter of the EMFN model has been established.


## 2.4. A Modified Method for Complex Fracture Network

The functionkTm,L,lmin describes the relationship between macroscopic hydraulic parameters and the microscopic parameters in EMFN. It is obtained from data fitting for the mesh that does not include any fracture, which means that the nodes used to generate the Delaunay-based mesh satisfy the “random node saturation procedure” in any part of the research domain. In other words, the density of chaos-oriented interfaces is homogeneous, and the influence of a slight disturbance in node density is negligible. However, when it comes to complex fracture networks (such as those shown in Figure 8, where the green disks are fractures and the blue cube is the research domain), the effect of additional refinement nodes cannot be ignored because they are too many. Where the refinement nodes are added, the density of the interface will be higher than that of other areas where they do not contain refinement nodes. Refinement nodes exist only in the fracture plane. Under this feature, according to the formation rule of Delaunay’s tetrahedron and the EMFN model, the equivalent permeability of the matrix will be affected only within the distance of one tetrahedron around the refinement node, and the transmissivity within the fracture plane will not be affected. The denser the interface network is, the higher the matrix network equivalent permeability is.Figure 8
Illustration for a complex fracture network (Nfr=10, 40, 80).
(a)
(b)
(c)In light of the problem, a practical strategy to reduce the influence of refinement nodes on the equivalent permeability of the matrix network is developed. We mark all the refinement nodes added by Step 3 introduced in Section2.2 and then find out in all the matrix interfaces that at least one of their vertex nodes is the refinement node. Finally, the conductivity of these interfaces is reduced because they are denser than we expected. Through a series of tests, a pragmatic reduction coefficient, 0.55, was determined. That is, if Tm is the transmissivity of matrix interfaces used in (14), then the transmissivity of these dense matrix interfaces is 0.45 Tm. We cannot prove mathematically that the reduction of transmissivity is correct, but it works.A series of experiments were carried out to verify the effectiveness of the reduction strategy. We designed a research domain thatlmin/L=0.02, and in which the centers of fractures are randomly distributed, and their orientation is random and isotropic with periodic boundary conditions (i.e., the example in Figure 8). The radius of fractures is identical and equal to 0.25 L. Under these conditions, samples with different fracture densities were tested. Fracture densities are assessed using the number Nfr of fractures inside the research domain. In this test, all the transmissivities of interfaces are identical and equal to Tm. Thus, although fractures are indeed described in the mesh, the fracture interfaces are no different from the matrix interfaces. So the influence of nonhomogeneous interface density and the effect of transmissivity reduction strategy can be easily known through steady-state flow simulation. We define the hydraulic conductivity of the matrix calculated by (14) as K0 and the hydraulic conductivity of fractured media obtained through EMFN as K. As shown in Figure 9, the black squares are K/K0 under the reduction strategy, which does not change significantly with the increase in fracture number (more specifically, the rise of refinement node). However, the results (red spots) obtained from the samples that not have a reduction procedure are obviously larger. The percentages of refinement nodes in total nodes are also noted in Figure 9.Figure 9
Effects of reduction strategy for samples with differentNfr.


## 2.5. Verification of the EMFN Model

In the EMFN model, either fracture interfaces or the equivalent matrix interfaces are all treated in the same way of the equivalent-fracture network, in which the “fracture interfaces” and the “matrix interfaces” are givenT independently; for fracture it is Tf, and for matrix it is Tm. Unfortunately, for flow problems in the discrete fracture porous media, it is hard to find an analytic solution of conductivity even for very simple conditions [27]. However, we can use the superposition principle in conductivity estimation to estimate the conductivity in some well-designed conditions. Two examples are given in this section. The first one is to test the conductivity of the rock a single rectangular fracture passing through the domain; the other one is to examine the conductivity tensor of the rock which has a disk fracture embedded with varying orientations in the centroid of the domain.When a rectangular fracture is parallel to the hydraulic gradientJ (sketching is depicted in Figure 10), its contribution to flow rate in the direction of J is easy to get by
(15)Qf=TflfJ,where Qf is the flow rate through fracture and lf is the trace length of the fracture. In this example, lf=60m, Tf=1×10−8m2/s, Tm=1×10−12m2/s, and L/lmin=50. Substituting (15) into (9) results in the x-direction as θfx=6×103. The transmission capacity of the matrix can be calculated by (9), (12), and (14) as θmx=2.18×102. The result of simulation using the EMFN model is θx=6.217×103.Figure 10
Single rectangular crack throughout the domain.
(a)
(b)The second example is illustrated in Figure11(a): a disk fracture (30 m in diameter) is embedded in the cubic domain, Tf=1×10−8m2/s, Tm=1×10−12m2/s, and L/lmin=50. The boundary condition is depicted in Figure 11(c), and a pair of boundaries perpendicular to the x-axis are given, respectively, constant hydraulic head h1=100m and h2=0m; others are corresponding linear variable hydraulic head boundaries. The center of the disk fracture coincides with the centroid of the cubic domain; the normal vector n of the fracture is always perpendicular to the z-axis. α is defined as the angle formed by the x-axis and n. The conductivity tensors when α equal to 0°, 15°, 30°, 45°, 60°, 75°, and 90° are calculated.Figure 11
Single circular fracture rotating in the center.
(a)
(b)
(c)The conductivity tensorK is the function of α and can be expressed as
(16)Kα=KxxαKxyα0KyxαKyyα000Kzzα.Whenα=0°,
(17)Kα=K1000K2000K3,where K1,K2,andK3 are, respectively, the coefficient of conductivity in the x,y,andz directions. Components in (16) can be written as
(18)Kxxα=K1cos2α+K2sin2α,Kyyα=K1sin2α+K2cos2α,Kxyα=Kyxα=K1−K2sinαcosα,Kzzα=K3.SubstitutingK1,K2,andK3 adopted from the numerical simulation when α=0° (here K1=2.111×10−12m/s, K2=K3=2.144×10−12m/s) into (16), (17), and (18), the relationships between α and components of K are received. K0 is defined as the coefficient of conductivity simulated in the same conditions but without any fracture (in this case, K0=2.111×10−12m/s and equal to K1, which means the fracture has no improvement on the conductivity in the direction parallel to n). The comparison between the numerical result using the EMFN model and the analytical result calculated by (16), (17), and (18) is plotted in Figure 12. One can see that the numerical results obtained are very close to the theoretical prediction values, which proves the validity of the EMFN method proposed in this paper.Figure 12
The comparison theoretical (solid line) and the numerical (symbol) result of components ofK.


## 3. Effects of Fracture Geometry on the Effective Conductivity

Different from general models based on the DFN model, the EMFN model is very suitable for studying the effects of the sparse fracture network on the effective conductivity of rock masses in which the contribution of the matrix to the macroscopic conductivity cannot be neglected. In this section, the effect of the hydraulic aperture, the relative position between fractures, and the connectivity of the fracture network on the conductivity of rock mass are studied by the EMFN model.


### 3.1. Hydraulic Aperture of Fracture

In this example, a single disk fracture is embedded in the cubic domain (L=100m). The fracture is parallel to the xoy-plane and has the same centroid as the cube (illustrated in Figure 13). The diameter of the disk is set as, respectively, 40 m, 50 m, and 60 m. Mesh L/lmin=50. Kx is defined as the coefficient of conductivity in the x direction, and K0 is the same coefficient but in the condition that no fracture was embedded. In this paper, the fracture is simplified to the disk plane. According to the cubic law, the aperture of fracture is proportional to the transmissivity in the plane:
(19)Tf=gbf312υ,where υ and g are, respectively, the coefficient of kinematic viscosity of water and the gravity acceleration. A series of simulations are conducted by varying the fracture transmissivity Tf and the diameter of fracture. Results are depicted in Figure 14, where the ordinate represents the relative hydraulic conductivity Kx/K0 and the abscissa is the ratio of fracture transmissivity and the transmissivity of the equivalent matrix-fracture network Tf/Tm. It can be observed that hydraulic conductivity increases until Tf/Tm grows to 105, after which effective conductivity in the x direction will remain constant and the increase in Tf has no more impact on Kx.Figure 13
Single horizontal circular fracture.Figure 14
Relationship betweenK/K0 and Tf/Tm.


### 3.2. Interaction between Two Fractures

From Section2.5, we know that a single fracture has almost no contribution to the macroscopic conductivity of rock in the same direction of its normal vector. However, the interaction of multiple fractures may have a complex impact on seepage behaviors of rock. Illustrated in Figure 15, two disk fractures (D=60m, Tf/Tm=104) are embedded in the research domain. For the mesh, L=100m and L/lmin=50. Fracture 1 is put at the center of the cubic domain, and its normal vector n1 is parallel to the z-axis; for fracture 2, the normal vector n2 is parallel to the x-axis, and the center of the fracture moves along the centerline of the cube in the x-direction. Thus, we can use the value of the x-coordinate (0m<x<100m) of the centroid node of the fracture to describe its spatial position. We propose to test the coefficient of conductivity in three directions (parallel to three axes), which are Kx,Ky,andKz. The change in relative conductivity in three directions is shown in Figure 16. One can observe that as fracture 2 moves, trends of Kx/K0, Ky/K0, and Kz/K0 are completely different.Figure 15
Illustration of the position of fractures and the contour plane.Figure 16
Relative permeabilityK/K0 in three directions.From the test in Section3.1, we obtained that Kx/K0=1.16 when there is only one single fracture (D=60m, Tf/Tm=104, and L/lmin=50). In this section, Kx/K0 also tends to 1.16 in the x direction when fracture 2 is moving away from fracture 1 and, notably, when fracture 2 is at the center of the cubic domain. Figure 17 illustrates the contours of the hydraulic head at the intersecting surface y=50m when the direction of J is parallel to the x-axis, which corresponds to the black squares in Figure 16. Obviously, the hydraulic head field in the cubic domain is symmetrical about this section. One can observe that the hydraulic head field is significantly influenced by the appearance of fracture 2, which makes the difference of Kx/K0. In fact, it can be seen from the distribution of the hydraulic head field on the fractures that although the individual fracture 2 does not contribute to the conductivity of rock mass in the x direction, it helps fracture 1 to obtain more flow recharge, making Kx/K0 much larger than the one when only fracture 1 is embedded. In particular, the hydraulic head field in Figures 17(a) and 17(d) are almost the same. In Figure 17(d), there is almost no hydraulic head change in fracture 2, so it barely helps the fluid flow in fracture 1, which explains the same value of Kx/K0 in Figure 16 at x=0m and x=50m. Moreover, in Figure 17, we can see that the change of the hydraulic head in fractures is much slower than in the rock matrix; in other words, it consumes much less energy of fluid flowing through the fracture than the rock matrix in the same distance, which is reasonable.Figure 17
Contour of the hydraulic head field aty=50m and on fracture plane.
(a)
(b)
(c)
(d)


### 3.3. Rock Mass with Multiple Fracture Clusters

Illustrated in Figure18 is a cubic domain with 3 fracture clusters embedded, and the distribution parameters are listed in Table 2. Following Baecher et al.’s model [28], fractures are supposed to be disk-shaped. The orientation of the normal vector of fractures obeys the Fisher distribution in each cluster, where κ is the dispersion parameter about the mean direction of fracture normal vectors [29]. The centroid position and size of fractures follow, respectively, uniform distribution and negative exponential distribution. Here, L=100m, L/lmin=50, and Tf/Tm=104; the hydraulic head upstream and downstream boundaries are, respectively, 100 m and 0 m.Figure 18
Contour of the hydraulic head field on fracture planes and matrix.
(a)
(b)
(c)
(d)Table 2
Parameters of fracture clusters.
Set
Amount of fractures
Normal direction
Mean diameter (m)
Dip direction/angle (deg)
κ
1
68
30/75
150
20
2
64
95/60
150
25
3
59
240/10
150
30To show the hydraulic head field clearly, the equivalent matrix-fractures are hidden. Figure18(a) is the hydraulic head in fractures. Obviously, the fracture network connects the upstream and downstream boundaries. Figure 18(c) is the hydraulic head contour in the rock matrix; the black thick solid line is the intersecting line formed by the contour plane and macroscopic fractures. To compare with the previous example, the geometric structure of the fracture network in Figures 18(b) and 18(d) is identical to Figures 18(a) and 18(c), but the fractures intersecting with plane x=50m are all eliminated, which means that the fracture network does not penetrate the domain in the x-direction. One can see that the connectivity of the fracture network significantly influences the hydraulic head field in fractured porous media. When fractures cut between each other and form a good percolation pathway, the distribution of the hydraulic gradient is relatively large along the path of fractures (the left side in Figure 18). However, for the nonpercolating situation (see Figure 18(b), the fracture did not cross the research domain in the x-direction); the change of the hydraulic head occurs mainly in the matrix (see Figure 18(d)). The relative conductivity Kx/K0 under different Tf/Tm is depicted in Figure 19 (the data of dark solid lines is taken from the example shown in Figure 18, and other lighter ones are obtained from ten parallel tests under the same fracture distributions). It can be seen from Figure 19 that the relative conductivity Kx/K0 of fractured rock mass increases with the increase in Tf/Tm. The difference is when Tf/Tm>104; Kx/K0 still increases and has a linear relationship with Tf/Tm in the percolating one, but Kx/K0 becomes almost a constant in the nonpercolating one. Corresponding to the situation of single fracture in Section 3.1, the fracture network that does not penetrate the entire porous media which can enhance the conductivity of the whole domain, but the conductivity of the rock matrix limits its upper limit.Figure 19
Change ofKx/K0 under different Tf/Tm of 10 specimens (dark solid lines are taken from the samples in Figure 18).


## 3.1. Hydraulic Aperture of Fracture

In this example, a single disk fracture is embedded in the cubic domain (L=100m). The fracture is parallel to the xoy-plane and has the same centroid as the cube (illustrated in Figure 13). The diameter of the disk is set as, respectively, 40 m, 50 m, and 60 m. Mesh L/lmin=50. Kx is defined as the coefficient of conductivity in the x direction, and K0 is the same coefficient but in the condition that no fracture was embedded. In this paper, the fracture is simplified to the disk plane. According to the cubic law, the aperture of fracture is proportional to the transmissivity in the plane:
(19)Tf=gbf312υ,where υ and g are, respectively, the coefficient of kinematic viscosity of water and the gravity acceleration. A series of simulations are conducted by varying the fracture transmissivity Tf and the diameter of fracture. Results are depicted in Figure 14, where the ordinate represents the relative hydraulic conductivity Kx/K0 and the abscissa is the ratio of fracture transmissivity and the transmissivity of the equivalent matrix-fracture network Tf/Tm. It can be observed that hydraulic conductivity increases until Tf/Tm grows to 105, after which effective conductivity in the x direction will remain constant and the increase in Tf has no more impact on Kx.Figure 13
Single horizontal circular fracture.Figure 14
Relationship betweenK/K0 and Tf/Tm.


## 3.2. Interaction between Two Fractures

From Section2.5, we know that a single fracture has almost no contribution to the macroscopic conductivity of rock in the same direction of its normal vector. However, the interaction of multiple fractures may have a complex impact on seepage behaviors of rock. Illustrated in Figure 15, two disk fractures (D=60m, Tf/Tm=104) are embedded in the research domain. For the mesh, L=100m and L/lmin=50. Fracture 1 is put at the center of the cubic domain, and its normal vector n1 is parallel to the z-axis; for fracture 2, the normal vector n2 is parallel to the x-axis, and the center of the fracture moves along the centerline of the cube in the x-direction. Thus, we can use the value of the x-coordinate (0m<x<100m) of the centroid node of the fracture to describe its spatial position. We propose to test the coefficient of conductivity in three directions (parallel to three axes), which are Kx,Ky,andKz. The change in relative conductivity in three directions is shown in Figure 16. One can observe that as fracture 2 moves, trends of Kx/K0, Ky/K0, and Kz/K0 are completely different.Figure 15
Illustration of the position of fractures and the contour plane.Figure 16
Relative permeabilityK/K0 in three directions.From the test in Section3.1, we obtained that Kx/K0=1.16 when there is only one single fracture (D=60m, Tf/Tm=104, and L/lmin=50). In this section, Kx/K0 also tends to 1.16 in the x direction when fracture 2 is moving away from fracture 1 and, notably, when fracture 2 is at the center of the cubic domain. Figure 17 illustrates the contours of the hydraulic head at the intersecting surface y=50m when the direction of J is parallel to the x-axis, which corresponds to the black squares in Figure 16. Obviously, the hydraulic head field in the cubic domain is symmetrical about this section. One can observe that the hydraulic head field is significantly influenced by the appearance of fracture 2, which makes the difference of Kx/K0. In fact, it can be seen from the distribution of the hydraulic head field on the fractures that although the individual fracture 2 does not contribute to the conductivity of rock mass in the x direction, it helps fracture 1 to obtain more flow recharge, making Kx/K0 much larger than the one when only fracture 1 is embedded. In particular, the hydraulic head field in Figures 17(a) and 17(d) are almost the same. In Figure 17(d), there is almost no hydraulic head change in fracture 2, so it barely helps the fluid flow in fracture 1, which explains the same value of Kx/K0 in Figure 16 at x=0m and x=50m. Moreover, in Figure 17, we can see that the change of the hydraulic head in fractures is much slower than in the rock matrix; in other words, it consumes much less energy of fluid flowing through the fracture than the rock matrix in the same distance, which is reasonable.Figure 17
Contour of the hydraulic head field aty=50m and on fracture plane.
(a)
(b)
(c)
(d)


## 3.3. Rock Mass with Multiple Fracture Clusters

Illustrated in Figure18 is a cubic domain with 3 fracture clusters embedded, and the distribution parameters are listed in Table 2. Following Baecher et al.’s model [28], fractures are supposed to be disk-shaped. The orientation of the normal vector of fractures obeys the Fisher distribution in each cluster, where κ is the dispersion parameter about the mean direction of fracture normal vectors [29]. The centroid position and size of fractures follow, respectively, uniform distribution and negative exponential distribution. Here, L=100m, L/lmin=50, and Tf/Tm=104; the hydraulic head upstream and downstream boundaries are, respectively, 100 m and 0 m.Figure 18
Contour of the hydraulic head field on fracture planes and matrix.
(a)
(b)
(c)
(d)Table 2
Parameters of fracture clusters.
Set
Amount of fractures
Normal direction
Mean diameter (m)
Dip direction/angle (deg)
κ
1
68
30/75
150
20
2
64
95/60
150
25
3
59
240/10
150
30To show the hydraulic head field clearly, the equivalent matrix-fractures are hidden. Figure18(a) is the hydraulic head in fractures. Obviously, the fracture network connects the upstream and downstream boundaries. Figure 18(c) is the hydraulic head contour in the rock matrix; the black thick solid line is the intersecting line formed by the contour plane and macroscopic fractures. To compare with the previous example, the geometric structure of the fracture network in Figures 18(b) and 18(d) is identical to Figures 18(a) and 18(c), but the fractures intersecting with plane x=50m are all eliminated, which means that the fracture network does not penetrate the domain in the x-direction. One can see that the connectivity of the fracture network significantly influences the hydraulic head field in fractured porous media. When fractures cut between each other and form a good percolation pathway, the distribution of the hydraulic gradient is relatively large along the path of fractures (the left side in Figure 18). However, for the nonpercolating situation (see Figure 18(b), the fracture did not cross the research domain in the x-direction); the change of the hydraulic head occurs mainly in the matrix (see Figure 18(d)). The relative conductivity Kx/K0 under different Tf/Tm is depicted in Figure 19 (the data of dark solid lines is taken from the example shown in Figure 18, and other lighter ones are obtained from ten parallel tests under the same fracture distributions). It can be seen from Figure 19 that the relative conductivity Kx/K0 of fractured rock mass increases with the increase in Tf/Tm. The difference is when Tf/Tm>104; Kx/K0 still increases and has a linear relationship with Tf/Tm in the percolating one, but Kx/K0 becomes almost a constant in the nonpercolating one. Corresponding to the situation of single fracture in Section 3.1, the fracture network that does not penetrate the entire porous media which can enhance the conductivity of the whole domain, but the conductivity of the rock matrix limits its upper limit.Figure 19
Change ofKx/K0 under different Tf/Tm of 10 specimens (dark solid lines are taken from the samples in Figure 18).


## 4. Conclusion and Perspective

In this study, an original method for simulating fluid flow in three-dimensional fractured porous media is proposed. The porous matrix and fractures are hydraulically replaced by an equivalent matrix-fracture network (EMFN). The finite element method is used to solve steady-state seepage flow in discretized triangular mesh. This method can naturally describe geometrical properties of each fracture and the contribution on equivalent conductivity of the fracture network and porous matrix.In the proposed method, two major obstacles have been overcome: one is the 3D mesh generation techniques based on the constrained Delaunay algorithm that can describe the geometry and topology of the fracture network and the rock matrix surrounding it; the other one is the relationship between macroscopic hydraulic parameters and the microscopic parameters in EMFN which ensures that the equivalent conductivity of fractured porous media is correctly simulated. The validity of the original method is verified through two simple examples from the perspective of equivalent conductivity and conductivity tensor. Moreover, the code is subsequently applied to study the effects of a single fracture aperture and the relative position of two fractures on the equivalent conductivity of fractured porous media. It is found that isolated fracture and fracture clusters surrounded by the porous media can also significantly influence the hydraulic head field distribution in its vicinity and increase the equivalent conductivity as those percolating ones. This effect depends on the properties of fractures individually and their relative position, which means the effective conductivity of fractured porous media cannot be obtained by simply accumulating the contributions from the porous matrix and each fracture. Finally, the code is used to simulate the steady-state flow in rock with 3 fracture clusters, and one controlled test is constructed. The results show that for fractured porous media with similar geometrical statistical parameters, the contribution of the matrix to equivalent conductivity becomes smaller as the conductivity of fractures increases. But for those where the fracture network does not penetrate the entire media, the equivalent conductivity will eventually reach an upper limit and will not increase as the conductivity of the fracture network increases.Up to now, the original method has been successfully used in steady-state flow simulation in fractured porous media. However, there are still several areas that need to be improved in the future. The density of the mesh is required to be uniformly distributed in space in the current method, which is computationally demanding when the fracture network is sparse and could cause the higher evaluation of matrix conductivity when the fracture network is too dense. Therefore, it is important to develop a high-quality variable density mesh generation technology. Moreover, the method can be further developed for hydraulic-mechanical coupling problems based on 3D DFN representation (as is done in the previous work for the 2D problem [22]), and for the stability analyses for the fractured rock [30–32].



---

*Source: 9736729-2019-01-03.xml*