# Hybrid Analytical and Numerical Approach for Modeling Fluid Flow in Simplified Three-Dimensional Fracture Networks

Authors: D. Roubinet; S. Demirel; E. B. Voytek; X. Wang; J. Irving
Journal: Geofluids
 (2020)

---

## Abstract

Modeling fluid flow in three-dimensional fracture networks is required in a wide variety of applications related to fractured rocks. Numerical approaches developed for this purpose rely on either simplified representations of the physics of the considered problem using mesh-free methods at the fracture scale or complex meshing of the studied systems resulting in considerable computational costs. Here, we derive an alternative approach that does not rely on a full meshing of the fracture network yet maintains an accurate representation of the modeled physical processes. This is done by considering simplified fracture networks in which the fractures are represented as rectangles that are divided into rectangular subfractures such that the fracture intersections are defined on the borders of these subfractures. Two-dimensional analytical solutions for the Darcy-scale flow problem are utilized at the subfracture scale and coupled at the fracture-network scale through discretization nodes located on the subfracture borders. We investigate the impact of parameters related to the location and number of the discretization nodes on the results obtained, and we compare our results with those calculated using reference solutions, which are an analytical solution for simple configurations and a standard finite-element modeling approach for complex configurations. This work represents a first step towards the development of 3D hybrid analytical and numerical approaches where the impact of the surrounding matrix will be eventually considered.


---

## Body

## 1. Introduction

Modeling fluid flow in the subsurface is required for numerous research fields and applications (e.g., [1–5]). This critical task is especially challenging in fractured rocks, as a large amount of work has shown that standard continuum representations cannot be used for most of the fractured porous domains encountered in the natural environment (e.g., [6], and references therein). This has led to the development of alternative subsurface representations where the interconnected fractures are represented individually as one- or two-dimensional discrete elements in two- and three-dimensional domains (e.g., [7–9]). These representations are well suited to capturing the impact of fracture-network properties on large-scale observations (e.g., [10–13]) and can be conditioned to data that are difficult to incorporate into continuum models (e.g., [14, 15]).Numerical modeling in fractured media using fracture-network representations is widespread and includes applications to multiphase reservoir flow (e.g., [16, 17]), coupled fluid flow, and solute transport for hydrogeological problems (e.g., [10, 13]), as well as electric current propagation in order to investigate fracture-network characteristics using geoelectrical methods (e.g., [18–20]). In this regard, when it is necessary to take into account fracture-matrix exchanges, such as when modeling electric current propagation or solute transport subject to matrix diffusion, several methods enable us to couple an explicit representation of the fractures with the surrounding matrix (e.g., [21–23]). In numerous hydrogeological studies, however, fluid flow is assumed to only occur in the interconnected fracture network, with the underlying assumption being that the surrounding matrix is impervious (e.g., [24–27]).Independent of whether fracture-matrix exchanges are ignored or considered, several strategies exist for simulating fluid flow in 3D fracture networks, which rely on different assumptions about the modeled processes at the fracture scale and the geometrical complexity of the fractures and fracture networks. The most cost-effective numerical approaches consider 1D fluid flow at the fracture scale, where the fractures are represented as pipes or rectangles, and the fracture-network discretization only requires nodes at the fracture extremities and intersections when fracture-matrix exchanges are not considered (e.g., [8, 28]). Considering rectangular fractures, the latter approach has been extended to 2D fluid flow at the fracture scale by discretizing each fracture with structured quadrilateral meshes and focusing on simple fracture-network configurations and intersections (e.g., [17, 29]). Finally, numerical methods have been developed for complex fracture networks relying upon more advanced representations of (i) the fracture shape, commonly using ellipses or circles, and (ii) the fracture-network geometry, by having no restriction concerning the corresponding fracture intersections (e.g., [24, 30]). These improvements were made at the expense of the computational cost of the corresponding simulations, which require the use of unstructured meshes associated with sophisticated meshing algorithms (e.g., [31]).In this paper, we investigate how analytical solutions at the fracture scale can be used to simulate Darcy-scale fluid flow processes in 3D fracture networks. To this end, each fracture is represented as a rectangle and divided into subfractures that are defined by the fracture intersections at the fracture-network scale. 2D analytical solutions are used in the subfractures and coupled through discretization nodes located on the subfracture borders, by enforcing continuity and mass conservation laws for the hydraulic head and fluid flow, respectively. The surface area of the considered 2D rectangular elements is not meshed, while keeping a 2D representation of the fluid flow in each element.After presenting this new modeling approach with two different means of defining the subfractures, we explain some restrictions related to the geometrical complexity of the considered fracture networks, and we provide the hydraulic head distribution at the fracture scale and the flow rate exiting the domain at the fracture-network scale for simple and complex fracture networks, respectively. In the latter case, the impact of the model parameters on the results is analyzed and these results are compared to those obtained using a standard finite-element approach.


## 2. Methodology

### 2.1. Governing Equations

Modeling fluid flow in 3D fracture networks is generally done by representing the fractures as 2D elements in which the flow rate is averaged over the fracture apertureb. When considering a Newtonian fluid at low Reynolds number, the flow within the fracture is governed by the Stokes equations and the corresponding average flow rate is expressed with the Darcy equation. In the context of an incompressible fluid under steady-state conditions, this results in the following expressions for mass conservation and the average flow rate per unit length J for each fracture (e.g., [32, 33]):
(1)∇·J=0,J=−T∇h.In equations (1), the hydraulic fracture transmissivity is given by T=ρgb3/12μ when assuming local Poiseuille law, with ρ and μ the fluid density and dynamic viscosity, respectively, g the acceleration of gravity, and h the hydraulic head.


### 2.2. Modeling Strategy

Equations (1) must be verified for each fracture of the considered fracture network. When considering a 2D representation of the fluid flow in the fractures, this is usually done by expressing these equations in the local coordinates of the fracture and solving them numerically by discretizing the fracture (see Introduction for references). In this work, we would like to avoid such a discretization by using analytical solutions of the mass conservation equation (1). To do so, we represent the fractures as 2D rectangles that are divided into rectangular subfractures such that the intersection lines between fractures are defined on the borders of these subfractures. This enables us to avoid the presence of intersections inside the 2D rectangular elements for which the analytical solutions are defined. For a simple example of two intersecting fractures, Figure 1 illustrates the corresponding subfractures that are determined with two different definitions. In the first case, the subfractures are defined such that the intersection line between the two main initial fractures corresponds exactly with the border of at least one subfracture. In the second case, the subfractures are defined such that the intersection line between the two main initial fractures is on the border of at least one subfracture but does not need to be fully defined on this border. We name these two definitions Inter1 and Inter2, respectively, where the latter allows us to reduce the number of subfractures compared with the former. Note that both definitions require, as mentioned previously, that (i) each fracture be represented by a 2D rectangle and (ii) the intersection lines between fractures have the same orientation as the fracture borders.Figure 1
Example of (a) two intersecting fractures that have been divided into (b) 8 subfractures with the methodInter1 and (c) 4 subfractures with the method Inter2. The intersection line between the two main initial fractures is shown in red.
(a)
(b)
(c)


### 2.3. Analytical Solution for the Subfractures

In each subfracture, the analytical solution of the mass conservation equation (1) is determined by considering a rectangle having length L, width H, and constant transmissivity T. In the domain Ω=xf,yf:0≤xf≤L,0≤yf≤H where xf and yf are the local Cartesian coordinates defined along the subfracture borders, the corresponding boundary-value problem (BVP) is expressed as
(2a)∂xf2hxf,yf+∂yf2hxf,yf=0,(2b)hxf,yf=gixf,yf,xf,yf∈Γi,where Γii=1,2,3,4 are the borders of the considered subfracture as shown in Figure 2.Figure 2
Illustration of a subfracture having lengthL, width H, local coordinates xf and yf, and discretized borders Γii=1,2,3,4 where the discretization nodes are represented with crosses.As described in AppendixA, the solution hx,y of BVP (2) is expressed as
(3a)hxf,yf=∑i=142li∑k=1∞sinkπχili∫0lisinkπziligizidzisinhkπχ~i/lisinhkl~iπ/li,with
(3b)li=H,i=1,2,L,i=3,4,l~i=L,i=1,2,H,i=3,4,χi=yf,i=1,2,xf,i=3,4,χ~i=L−xf,i=1,xf,i=2,H−yf,i=3,yf,i=4,zi=yf,i=1,2,xf,i=3,4.In expression (3), the functions gi describe the unknown hydraulic head h along the borders Γi, which are discretized as shown in Figure 2. Discretizing each border Γi into Ni nodes results in discretizing the corresponding unknown function gi into Ni discretized values denoted gipp=1,⋯,Ni. These values are used to enforce the flow conditions described next in Section 2.4, and they are the unknowns of the linear system defined in Section 2.5.


### 2.4. Flow Conditions for the Fracture Network

Solving equation (1) in 3D fracture networks is accomplished by enforcing the following flow conditions: (1) flow continuity at the interfaces between subfractures; (2) no-flow conditions on the fracture borders, as we consider here a surrounding matrix impervious to fluid flow; and (3) flow conservation at the fracture intersections. These three conditions are applied along the subfracture borders that are discretized as shown in Figure 2, implying that they can be mathematically expressed at each node m as follows:
(4)∑k∈IFmJkxkmΔm=0.In expression (4), IFm represents the set of subfractures that share node m; Jkxkm is the average flow rate per unit length leaving fracture k at position xkm=xkm,ykm, i.e., the position of node m in fracture k; and Δm is the length associated with the discretized node m. This parameter is defined as Δm=lfm/Nim with lfm the length of the fracture on which node m is located and Nim the number of discretized nodes on this fracture, which is determined from the following expression: Nim=min1,intlfm/Δdisc. In this expression, the function minx,y returns the minimum value between x and y and the function intx rounds x down to the integer below x. From this definition, there is at least one discretized node on each fracture segment and the minimum length associated with each discretization node is the length Δdisc.


### 2.5. Final Linear System

Expressions for the average flow rates required in (4) are given in Appendix B for each border of a given subfracture. These expressions depend linearly on the unknown discretized values of hydraulic head gip along each subfracture border. Writing the flow condition (4) for each discretization node leads to a linear system Ax=b where x is the vector of unknown values gip. After solving this linear system by considering Dirichlet or Neumann boundary conditions on the domain borders, the distributions of hydraulic head and average flow rate in the fracture network are defined using expression (3) and the Darcy equation (1), respectively, for each subfracture.


## 2.1. Governing Equations

Modeling fluid flow in 3D fracture networks is generally done by representing the fractures as 2D elements in which the flow rate is averaged over the fracture apertureb. When considering a Newtonian fluid at low Reynolds number, the flow within the fracture is governed by the Stokes equations and the corresponding average flow rate is expressed with the Darcy equation. In the context of an incompressible fluid under steady-state conditions, this results in the following expressions for mass conservation and the average flow rate per unit length J for each fracture (e.g., [32, 33]):
(1)∇·J=0,J=−T∇h.In equations (1), the hydraulic fracture transmissivity is given by T=ρgb3/12μ when assuming local Poiseuille law, with ρ and μ the fluid density and dynamic viscosity, respectively, g the acceleration of gravity, and h the hydraulic head.


## 2.2. Modeling Strategy

Equations (1) must be verified for each fracture of the considered fracture network. When considering a 2D representation of the fluid flow in the fractures, this is usually done by expressing these equations in the local coordinates of the fracture and solving them numerically by discretizing the fracture (see Introduction for references). In this work, we would like to avoid such a discretization by using analytical solutions of the mass conservation equation (1). To do so, we represent the fractures as 2D rectangles that are divided into rectangular subfractures such that the intersection lines between fractures are defined on the borders of these subfractures. This enables us to avoid the presence of intersections inside the 2D rectangular elements for which the analytical solutions are defined. For a simple example of two intersecting fractures, Figure 1 illustrates the corresponding subfractures that are determined with two different definitions. In the first case, the subfractures are defined such that the intersection line between the two main initial fractures corresponds exactly with the border of at least one subfracture. In the second case, the subfractures are defined such that the intersection line between the two main initial fractures is on the border of at least one subfracture but does not need to be fully defined on this border. We name these two definitions Inter1 and Inter2, respectively, where the latter allows us to reduce the number of subfractures compared with the former. Note that both definitions require, as mentioned previously, that (i) each fracture be represented by a 2D rectangle and (ii) the intersection lines between fractures have the same orientation as the fracture borders.Figure 1
Example of (a) two intersecting fractures that have been divided into (b) 8 subfractures with the methodInter1 and (c) 4 subfractures with the method Inter2. The intersection line between the two main initial fractures is shown in red.
(a)
(b)
(c)


## 2.3. Analytical Solution for the Subfractures

In each subfracture, the analytical solution of the mass conservation equation (1) is determined by considering a rectangle having length L, width H, and constant transmissivity T. In the domain Ω=xf,yf:0≤xf≤L,0≤yf≤H where xf and yf are the local Cartesian coordinates defined along the subfracture borders, the corresponding boundary-value problem (BVP) is expressed as
(2a)∂xf2hxf,yf+∂yf2hxf,yf=0,(2b)hxf,yf=gixf,yf,xf,yf∈Γi,where Γii=1,2,3,4 are the borders of the considered subfracture as shown in Figure 2.Figure 2
Illustration of a subfracture having lengthL, width H, local coordinates xf and yf, and discretized borders Γii=1,2,3,4 where the discretization nodes are represented with crosses.As described in AppendixA, the solution hx,y of BVP (2) is expressed as
(3a)hxf,yf=∑i=142li∑k=1∞sinkπχili∫0lisinkπziligizidzisinhkπχ~i/lisinhkl~iπ/li,with
(3b)li=H,i=1,2,L,i=3,4,l~i=L,i=1,2,H,i=3,4,χi=yf,i=1,2,xf,i=3,4,χ~i=L−xf,i=1,xf,i=2,H−yf,i=3,yf,i=4,zi=yf,i=1,2,xf,i=3,4.In expression (3), the functions gi describe the unknown hydraulic head h along the borders Γi, which are discretized as shown in Figure 2. Discretizing each border Γi into Ni nodes results in discretizing the corresponding unknown function gi into Ni discretized values denoted gipp=1,⋯,Ni. These values are used to enforce the flow conditions described next in Section 2.4, and they are the unknowns of the linear system defined in Section 2.5.


## 2.4. Flow Conditions for the Fracture Network

Solving equation (1) in 3D fracture networks is accomplished by enforcing the following flow conditions: (1) flow continuity at the interfaces between subfractures; (2) no-flow conditions on the fracture borders, as we consider here a surrounding matrix impervious to fluid flow; and (3) flow conservation at the fracture intersections. These three conditions are applied along the subfracture borders that are discretized as shown in Figure 2, implying that they can be mathematically expressed at each node m as follows:
(4)∑k∈IFmJkxkmΔm=0.In expression (4), IFm represents the set of subfractures that share node m; Jkxkm is the average flow rate per unit length leaving fracture k at position xkm=xkm,ykm, i.e., the position of node m in fracture k; and Δm is the length associated with the discretized node m. This parameter is defined as Δm=lfm/Nim with lfm the length of the fracture on which node m is located and Nim the number of discretized nodes on this fracture, which is determined from the following expression: Nim=min1,intlfm/Δdisc. In this expression, the function minx,y returns the minimum value between x and y and the function intx rounds x down to the integer below x. From this definition, there is at least one discretized node on each fracture segment and the minimum length associated with each discretization node is the length Δdisc.


## 2.5. Final Linear System

Expressions for the average flow rates required in (4) are given in Appendix B for each border of a given subfracture. These expressions depend linearly on the unknown discretized values of hydraulic head gip along each subfracture border. Writing the flow condition (4) for each discretization node leads to a linear system Ax=b where x is the vector of unknown values gip. After solving this linear system by considering Dirichlet or Neumann boundary conditions on the domain borders, the distributions of hydraulic head and average flow rate in the fracture network are defined using expression (3) and the Darcy equation (1), respectively, for each subfracture.


## 3. Method Analysis and Applications

The modeling approach described in Section2 is tested for simple (Section 3.1) and complex (Section 3.2) fracture configurations. In Section 3.1, the study of the hydraulic head distribution in a single fracture for various values of the discretization parameters leads to define a new model parameter denoted ε, for which various definitions are tested on a two-intersecting fracture configuration. In Section 3.2, the impact of the discretization parameters is analyzed for complex fracture-network configurations, and the results obtained with the two subfracture definitions are compared with those coming from a reference solution.


### 3.1. Method Analysis with Simple Configurations

We consider the simple fracture configurations presented in Figure3 that correspond to a single fracture (Config_1a) and two intersecting fractures (Config_1b) in 10×10×10 m3 domains. The fracture aperture is set to 10−3 m and we assume Dirichlet hydraulic head boundary conditions equal to 1 m and 0 m on the left and right sides of the domain, respectively, and Dirichlet boundary conditions varying linearly between these two values along the top, bottom, back, and front sides. The hydraulic head distribution in the single fracture can be described with the analytical expression hx,y=1−x/L that we use with L=10 m as the reference solution for Config_1a, whereas we use a standard finite-element approach as the reference solution for Config_1b and the configurations considered next in Section 3.2.Figure 3
Simple fracture configurations considered.
(a)
Config_1a
(b)
Config_1bUsing the reference solutions described above, we define the relative error in the hydraulic headEh1x,y=∣hx,y−hrefx,y∣/hrefx,y×100 with hx,y and hrefx,y as the hydraulic heads at position x,y computed with our hybrid approach and the reference solution, respectively. We also define Eh2x as the average error at each position x of the relative errors Eh1 computed along position y. Figure 4 shows these errors along the relative distance xL=x/L for Config_1a with various values of the hybrid model parameter Δdisc going from 5 to 0.1 m. As explained in Section 2.4, this parameter defines the length Δm associated with each node m that is used in equation (4). Figure 4 shows that decreasing the discretization parameter Δdisc leads to decrease both errors Eh1 and Eh2 with values below 2% and 1% when Δdisc≤1 m and 0.5 m, respectively, except for some instabilities of Eh1 that are observed when xL is located at the fracture extremities (i.e., xL close to 0 and 1 in Figure 4(a)). Conversely, when Δdisc=5 m, large values of Eh1 and Eh2 are observed because only two discretization nodes are used to describe the hydraulic head distribution along the fracture borders.Figure 4
Errors in hydraulic head along the relative fracture lengthxL=x/L for the fracture configuration Config_1a (Figure 3(a)) with the discretization length Δdisc set to 5, 1, 0.5, and 0.1 m. The symbols represent (a) the relative error Eh1 at each position x,y and (b) the relative error Eh2 averaged over position y. The dashed and dotted horizontal black lines represent errors of 1 and 2%, respectively.
(a)
(b)The presence of instabilities of the hydraulic head shown at the fracture extremities in Figure4(a) leads to define the parameter ε, which corresponds to a small distance around the fracture extremities at which the averaged fluid flow rates are computed (Appendix B) and the fluid flow conservation is enforced (equation (4)). In order to evaluate the most appropriate definition of ε, we evaluate the relative error in the fluid flow rate that exits the right side of the domain for the configurations presented in Figure 3. This error is denoted Ex and is computed by considering as the reference solution for Config_1a the analytical expression of the flow rate Q=−TH∇h with T and H as the fracture transmissivity and width, respectively, and ∇h as the hydraulic head gradient between the left and right sides of the domain. As mentioned before, the reference solution for Config_1b is computed with a standard finite-element approach. Figure 4(a) shows that the location and amplitude of the instabilities of the hydraulic head depend on the fracture discretization such that the amplitude increases and the location is closer to the fracture extremities when the discretization parameter Δdisc decreases. In order to take into account this observation, we consider that the definition of ε should be expressed as a function of Δdisc and we choose the linear relationship ε=αΔdisc with α being an adjusting coefficient. This definition is tested by evaluating the impact of various values of α on the error in flow rate Ex for the configuration Config_1a. The results displayed in Figure 5(a) show that Ex quickly decreases to values smaller than 1% when increasing the characteristic discretization length l=L/Δdisc with α set to 0.5, 1, and 1.5. That being said, using α=1 results in a quicker and monotonic decrease in Ex, implying that we will use this value for the rest of this work.Figure 5
Error in the exiting flow rateEx (%) along the characteristic discretization length l=L/Δdisc. (a) The relationship ε=αΔdisc is tested for various values of the coefficient α with the fracture configuration Config_1a (Figure 3(a)). (b, c) Various definitions of ε are tested for the fracture configuration Config_1b (Figure 3(b)) using the methods (b) Inter1 and (c) Inter2 to define the subfractures of the system. The dashed and dotted horizontal black lines represent errors of 1 and 2%, respectively.
(a)
Config_1a
(b)
Config_1b, Inter1
(c)
Config_1b, Inter2In the previous configuration, the size of the single fracture and the values of the discretization parameterΔdisc are such that Δm=Δdisc because each border of the considered fracture could be exactly divided into segments of length Δdisc. However, for more complex configurations, the fracture borders might be divided into segments of slightly larger length than Δdisc. In this case, the parameter Δm that is used in equation (4) and the discretization parameter Δdisc that is used to define Δm might be different. In addition, Δm might not be constant at the fracture scale because the value of Δm could be different from a fracture border to the other. In order to evaluate the best definition of ε in this case, we compute the relative error Ex for the configuration Config_1b (Figure 3(b)) considering the following definitions of ε: (i) ε=Δdisc, (ii) ε=Δav_fractm, and (iii) ε=Δav_borderm, where Δav_fractm and Δav_borderm correspond to the value of Δm averaged at the fracture and border scale, respectively. Figures 5(b) and 5(c) show the results that are obtained using the methods Inter1 and Inter2 to define the subfractures of the system. These results show that averaging Δm at the border scale leads to a more accurate description of the exiting flow rate with errors smaller than 2 and 1% using the methods Inter1 and Inter2, respectively. Consequently, the definition of ε=Δav_borderm will be used in the next section of our study.


### 3.2. Application to Fracture Networks

We wish now to apply the modeling approach presented in Section2 and the parameter definitions provided in Section 3.1 to various fracture networks. These networks must respect the geometry restrictions explained in Section 2 implying that the intersection lines between fractures and subfractures must be parallel to the fracture borders. This condition is fulfilled by the deterministic and random fracture networks that are described below. Our hybrid analytical and numerical method is tested on these networks for various hydraulic conditions and fluid flow directions, and the corresponding results are shown below and discussed in Section 4.In the line of the deterministic layered fracture models presented in Doolin and Mauldon [34], we consider two sets of joint networks denoted Config_2 and Config_3 that are shown in Figure 6. These networks, consisting of strata-bound, bed-normal joint sets separated by bedding planes, are representative configurations of joint systems in layered carbonate rocks [35]. The size of the models is 2×2×0.1 m3. In Config_2, three configurations of a double-layered joint network with different spacing contrast between the two layers are examined (Figures 6(a)–6(c)). The joint spacing is 0.3 m in the top and bottom layers of Config_2a. In Config_2b and Config_2c, the joint spacing is 0.3 m for the top layer joint set, while equal to 0.6 m for the bottom layer joint set. The difference between the two models is the joint locations in the bottom layer. In Config_3, we test three double-layered joint models with different overlap patterns between joint sets in the top and bottom layers (Figures 6(d)–6(f)). The overlap between joint sets in the top and bottom layers is 0.6 m for Config_3a, 0.3 m for Config_3b, and 0 m for Config_3c. To calculate flow through the joint networks, the fracture aperture is set to 10−3 m and the hydraulic heads are set to 1 m and 0 m on the top and bottom sides of the domain, respectively, with no-flow conditions enforced along the other borders.Figure 6
Deterministic fracture networks representing strata-bound, bed-normal joint sets in layered carbonate rocks. Joint networks considering (a–c) spacing contrast in adjacent beds and (d–f) overlapping pattern of joint sets in adjacent beds.
(a)
Config_2a
(b)
Config_2b
(c)
Config_2c
(d)
Config_3a
(e)
Config_3b
(f)
Config_3cIn relation with previous studies describing fractured aquifers and reservoirs, we also consider the sets of random fracture networksConfig_4 and Config_5 shown in Figure 7 with the percolation parameter p=∑iNf0.5li3/L3, with Nf being the number of fractures, li the length of fracture i, and L the domain size [24, 30, 36, 37]. In these 10×10×10 m3 domains, the center of each fracture is determined by drawing the corresponding coordinates x, y, and z from a uniform distribution ranging from 0 to the domain size L=10 m. The fracture aperture is set to 10−3 m, and the hydraulic heads are set to 1 m and 0 m on the left and right sides of the domain, respectively, and vary linearly between these two values along the other borders. The fracture networks Config_4 (Figures 7(a)–7(c)) correspond to orthogonal square fractures of size 5 m with parameter p set to 0.2 (Config_4a), 0.25 (Config_4b), and 0.3 (Config_4c). For the fracture networks Config_5 (Figures 7(d) and 7(e)), we consider the following characteristics: (i) the fracture length is uniformly distributed from 0 to L, (ii) the fracture length is oriented along the x- or y-axis with equal probability, (iii) the ratio between the fracture length and width is set to 2, (iv) the angle between the fracture and the x-y plane is uniformly distributed between 0 and π, and (v) parameter p is set to 0.5 (Config_5a) and 0.55 (Config_5b). For both sets of networks Config_4 and Config_5, we only consider the fractures that are connected to the right side of the domain and we compute the flow rate exiting this side.Figure 7
Random fracture networks considering (a–c) orthogonal families of fractures and (d, e) random orientations for various values of the percolation parameterp.
(a)
Config_4a, p=0.2
(b)
Config_4b, p=0.25
(c)
Config_4c, p=0.3
(d)
Config_5a, p=0.5
(e)
Config_5b, p=0.55Our modeling approach is tested with these fracture networks by defining the relative errorsEz and Ex in the flow rate exiting the bottom and right side of the domain, respectively. As before, we use a standard finite-element method as the reference solution. Figure 8 shows Ez for the sets of fracture networks Config_2 and Config_3 and Ex for Config_4 and Config_5 along the characteristic discretization length l=L/Δdisc with L=2 and 10 m for Config_2-Config_3 and Config_4-Config_5, respectively. These results are analyzed and discussed in the following section.Figure 8
Error in the exiting flow rates (a, b)Ez (%) and (c–f) Ex (%) along the characteristic length l=L/Δdisc with (a, b) L=2 m and (c–f) L=10 m for the fracture networks (a) Config_2 (Figures 6(a)–6(c)), (b) Config_3 (Figures 6(d)–6(f)), (c, d) Config_4 (Figures 7(a)–7(c)), and (e, f) Config_5 (Figures 7(d) and 7(e)). The method Inter1 is used in (b, full lines), (c), and (e) to define the subfractures of the systems and the method Inter2 in (b, dotted lines), (d), and (f). The dashed, dashdot, and dotted black lines with no symbol represent errors of 1, 2, and 10%, respectively.
(a)
Config_2
(b)
Config_3
(c)
Config_4, Inter1
(d)
Config_4, Inter2
(e)
Config_5, Inter1
(f)
Config_5, Inter2


## 3.1. Method Analysis with Simple Configurations

We consider the simple fracture configurations presented in Figure3 that correspond to a single fracture (Config_1a) and two intersecting fractures (Config_1b) in 10×10×10 m3 domains. The fracture aperture is set to 10−3 m and we assume Dirichlet hydraulic head boundary conditions equal to 1 m and 0 m on the left and right sides of the domain, respectively, and Dirichlet boundary conditions varying linearly between these two values along the top, bottom, back, and front sides. The hydraulic head distribution in the single fracture can be described with the analytical expression hx,y=1−x/L that we use with L=10 m as the reference solution for Config_1a, whereas we use a standard finite-element approach as the reference solution for Config_1b and the configurations considered next in Section 3.2.Figure 3
Simple fracture configurations considered.
(a)
Config_1a
(b)
Config_1bUsing the reference solutions described above, we define the relative error in the hydraulic headEh1x,y=∣hx,y−hrefx,y∣/hrefx,y×100 with hx,y and hrefx,y as the hydraulic heads at position x,y computed with our hybrid approach and the reference solution, respectively. We also define Eh2x as the average error at each position x of the relative errors Eh1 computed along position y. Figure 4 shows these errors along the relative distance xL=x/L for Config_1a with various values of the hybrid model parameter Δdisc going from 5 to 0.1 m. As explained in Section 2.4, this parameter defines the length Δm associated with each node m that is used in equation (4). Figure 4 shows that decreasing the discretization parameter Δdisc leads to decrease both errors Eh1 and Eh2 with values below 2% and 1% when Δdisc≤1 m and 0.5 m, respectively, except for some instabilities of Eh1 that are observed when xL is located at the fracture extremities (i.e., xL close to 0 and 1 in Figure 4(a)). Conversely, when Δdisc=5 m, large values of Eh1 and Eh2 are observed because only two discretization nodes are used to describe the hydraulic head distribution along the fracture borders.Figure 4
Errors in hydraulic head along the relative fracture lengthxL=x/L for the fracture configuration Config_1a (Figure 3(a)) with the discretization length Δdisc set to 5, 1, 0.5, and 0.1 m. The symbols represent (a) the relative error Eh1 at each position x,y and (b) the relative error Eh2 averaged over position y. The dashed and dotted horizontal black lines represent errors of 1 and 2%, respectively.
(a)
(b)The presence of instabilities of the hydraulic head shown at the fracture extremities in Figure4(a) leads to define the parameter ε, which corresponds to a small distance around the fracture extremities at which the averaged fluid flow rates are computed (Appendix B) and the fluid flow conservation is enforced (equation (4)). In order to evaluate the most appropriate definition of ε, we evaluate the relative error in the fluid flow rate that exits the right side of the domain for the configurations presented in Figure 3. This error is denoted Ex and is computed by considering as the reference solution for Config_1a the analytical expression of the flow rate Q=−TH∇h with T and H as the fracture transmissivity and width, respectively, and ∇h as the hydraulic head gradient between the left and right sides of the domain. As mentioned before, the reference solution for Config_1b is computed with a standard finite-element approach. Figure 4(a) shows that the location and amplitude of the instabilities of the hydraulic head depend on the fracture discretization such that the amplitude increases and the location is closer to the fracture extremities when the discretization parameter Δdisc decreases. In order to take into account this observation, we consider that the definition of ε should be expressed as a function of Δdisc and we choose the linear relationship ε=αΔdisc with α being an adjusting coefficient. This definition is tested by evaluating the impact of various values of α on the error in flow rate Ex for the configuration Config_1a. The results displayed in Figure 5(a) show that Ex quickly decreases to values smaller than 1% when increasing the characteristic discretization length l=L/Δdisc with α set to 0.5, 1, and 1.5. That being said, using α=1 results in a quicker and monotonic decrease in Ex, implying that we will use this value for the rest of this work.Figure 5
Error in the exiting flow rateEx (%) along the characteristic discretization length l=L/Δdisc. (a) The relationship ε=αΔdisc is tested for various values of the coefficient α with the fracture configuration Config_1a (Figure 3(a)). (b, c) Various definitions of ε are tested for the fracture configuration Config_1b (Figure 3(b)) using the methods (b) Inter1 and (c) Inter2 to define the subfractures of the system. The dashed and dotted horizontal black lines represent errors of 1 and 2%, respectively.
(a)
Config_1a
(b)
Config_1b, Inter1
(c)
Config_1b, Inter2In the previous configuration, the size of the single fracture and the values of the discretization parameterΔdisc are such that Δm=Δdisc because each border of the considered fracture could be exactly divided into segments of length Δdisc. However, for more complex configurations, the fracture borders might be divided into segments of slightly larger length than Δdisc. In this case, the parameter Δm that is used in equation (4) and the discretization parameter Δdisc that is used to define Δm might be different. In addition, Δm might not be constant at the fracture scale because the value of Δm could be different from a fracture border to the other. In order to evaluate the best definition of ε in this case, we compute the relative error Ex for the configuration Config_1b (Figure 3(b)) considering the following definitions of ε: (i) ε=Δdisc, (ii) ε=Δav_fractm, and (iii) ε=Δav_borderm, where Δav_fractm and Δav_borderm correspond to the value of Δm averaged at the fracture and border scale, respectively. Figures 5(b) and 5(c) show the results that are obtained using the methods Inter1 and Inter2 to define the subfractures of the system. These results show that averaging Δm at the border scale leads to a more accurate description of the exiting flow rate with errors smaller than 2 and 1% using the methods Inter1 and Inter2, respectively. Consequently, the definition of ε=Δav_borderm will be used in the next section of our study.


## 3.2. Application to Fracture Networks

We wish now to apply the modeling approach presented in Section2 and the parameter definitions provided in Section 3.1 to various fracture networks. These networks must respect the geometry restrictions explained in Section 2 implying that the intersection lines between fractures and subfractures must be parallel to the fracture borders. This condition is fulfilled by the deterministic and random fracture networks that are described below. Our hybrid analytical and numerical method is tested on these networks for various hydraulic conditions and fluid flow directions, and the corresponding results are shown below and discussed in Section 4.In the line of the deterministic layered fracture models presented in Doolin and Mauldon [34], we consider two sets of joint networks denoted Config_2 and Config_3 that are shown in Figure 6. These networks, consisting of strata-bound, bed-normal joint sets separated by bedding planes, are representative configurations of joint systems in layered carbonate rocks [35]. The size of the models is 2×2×0.1 m3. In Config_2, three configurations of a double-layered joint network with different spacing contrast between the two layers are examined (Figures 6(a)–6(c)). The joint spacing is 0.3 m in the top and bottom layers of Config_2a. In Config_2b and Config_2c, the joint spacing is 0.3 m for the top layer joint set, while equal to 0.6 m for the bottom layer joint set. The difference between the two models is the joint locations in the bottom layer. In Config_3, we test three double-layered joint models with different overlap patterns between joint sets in the top and bottom layers (Figures 6(d)–6(f)). The overlap between joint sets in the top and bottom layers is 0.6 m for Config_3a, 0.3 m for Config_3b, and 0 m for Config_3c. To calculate flow through the joint networks, the fracture aperture is set to 10−3 m and the hydraulic heads are set to 1 m and 0 m on the top and bottom sides of the domain, respectively, with no-flow conditions enforced along the other borders.Figure 6
Deterministic fracture networks representing strata-bound, bed-normal joint sets in layered carbonate rocks. Joint networks considering (a–c) spacing contrast in adjacent beds and (d–f) overlapping pattern of joint sets in adjacent beds.
(a)
Config_2a
(b)
Config_2b
(c)
Config_2c
(d)
Config_3a
(e)
Config_3b
(f)
Config_3cIn relation with previous studies describing fractured aquifers and reservoirs, we also consider the sets of random fracture networksConfig_4 and Config_5 shown in Figure 7 with the percolation parameter p=∑iNf0.5li3/L3, with Nf being the number of fractures, li the length of fracture i, and L the domain size [24, 30, 36, 37]. In these 10×10×10 m3 domains, the center of each fracture is determined by drawing the corresponding coordinates x, y, and z from a uniform distribution ranging from 0 to the domain size L=10 m. The fracture aperture is set to 10−3 m, and the hydraulic heads are set to 1 m and 0 m on the left and right sides of the domain, respectively, and vary linearly between these two values along the other borders. The fracture networks Config_4 (Figures 7(a)–7(c)) correspond to orthogonal square fractures of size 5 m with parameter p set to 0.2 (Config_4a), 0.25 (Config_4b), and 0.3 (Config_4c). For the fracture networks Config_5 (Figures 7(d) and 7(e)), we consider the following characteristics: (i) the fracture length is uniformly distributed from 0 to L, (ii) the fracture length is oriented along the x- or y-axis with equal probability, (iii) the ratio between the fracture length and width is set to 2, (iv) the angle between the fracture and the x-y plane is uniformly distributed between 0 and π, and (v) parameter p is set to 0.5 (Config_5a) and 0.55 (Config_5b). For both sets of networks Config_4 and Config_5, we only consider the fractures that are connected to the right side of the domain and we compute the flow rate exiting this side.Figure 7
Random fracture networks considering (a–c) orthogonal families of fractures and (d, e) random orientations for various values of the percolation parameterp.
(a)
Config_4a, p=0.2
(b)
Config_4b, p=0.25
(c)
Config_4c, p=0.3
(d)
Config_5a, p=0.5
(e)
Config_5b, p=0.55Our modeling approach is tested with these fracture networks by defining the relative errorsEz and Ex in the flow rate exiting the bottom and right side of the domain, respectively. As before, we use a standard finite-element method as the reference solution. Figure 8 shows Ez for the sets of fracture networks Config_2 and Config_3 and Ex for Config_4 and Config_5 along the characteristic discretization length l=L/Δdisc with L=2 and 10 m for Config_2-Config_3 and Config_4-Config_5, respectively. These results are analyzed and discussed in the following section.Figure 8
Error in the exiting flow rates (a, b)Ez (%) and (c–f) Ex (%) along the characteristic length l=L/Δdisc with (a, b) L=2 m and (c–f) L=10 m for the fracture networks (a) Config_2 (Figures 6(a)–6(c)), (b) Config_3 (Figures 6(d)–6(f)), (c, d) Config_4 (Figures 7(a)–7(c)), and (e, f) Config_5 (Figures 7(d) and 7(e)). The method Inter1 is used in (b, full lines), (c), and (e) to define the subfractures of the systems and the method Inter2 in (b, dotted lines), (d), and (f). The dashed, dashdot, and dotted black lines with no symbol represent errors of 1, 2, and 10%, respectively.
(a)
Config_2
(b)
Config_3
(c)
Config_4, Inter1
(d)
Config_4, Inter2
(e)
Config_5, Inter1
(f)
Config_5, Inter2


## 4. Results and Discussion

Figure8 shows the errors in the main exiting flow rate obtained with the hybrid analytical and numerical approach described in Section 2 for the fracture networks presented in Figures 6 and 7. The results presented in Figures 8(a) and 8(b) are used to evaluate the applicability of our method for modeling vertical fluid flow in deterministic joint models and the results presented in Figures 8(c)–8(f) for horizontal fluid flow in random fracture networks. For all these cases, the two definitions of the subfractures previously denoted Inter1 and Inter2 are tested, except in Figure 8(a) for which the tested fracture networks result in identical networks of subfractures when applying either method Inter1 or Inter2.In the latter case, Figure8(a) shows that the error in flow rate Ez quickly decreases for the configurations Config_2 when increasing the characteristic discretization length l and reaches values below 2% when l≥l1 and around 1% when l>l1 with l1=40. This critical value, which is represented as a dashed vertical black line in Figure 8(a), is obtained when Δdisc is set to 0.05 m, which corresponds to the width of the top and bottom layers of Config_2. This value also corresponds to the length of the smallest subfracture borders of the considered systems with an example of these borders shown in red in Figure 6(a). Consequently, when Δdisc≥0.05 m (i.e., l≤l1), only one node is used on these smallest borders, and when Δdisc<0.05 m (i.e., l>l1), more than one node is used on these borders. Thus, the results in Figure 8(a) show that the exiting flow rate in Config_2 is well described when using one node (or more) to discretize the smallest subfracture borders.In order to evaluate the impact of reducing the fracture length and modifying the overlap between the fractures located on the top and bottom layers, we consider the joint modelsConfig_3 presented in Figures 6(d)–6(f). These models are studied with the definitions of the subfractures Inter1 and Inter2 for which the error in flow rate Ez is represented in Figure 8(b) with full and dashed lines, respectively. Although the smallest subfracture length is equal to 0.05 m as before, the results in Figure 8(b) show that more than one node is required to discretize these smallest subfracture borders since Ez is larger than 2% when l=l1. Consequently, in comparison with Config_2, the decrease in the length of the fractures in Config_3 results in increasing the complexity of the fluid flow such that more discretization nodes are required for representing this flow. Figure 8(b) also shows that the number of nodes must be increased from Config_3a to Config_3c since Ez is smaller than 2% with Inter1 when l≥100, 400, and 800 for Config_3a, Config_3b, and Config_3c, respectively. From Config_3a to Config_3c, reducing the overlap between the fractures located on the top layer and those located on the bottom layer impacts the behavior of the fluid circulating from the top to the bottom fractures through the bedding plane and results in increasing the complexity of this flow. When the fractures fully overlap (Config_3a), the flow in the bedding plane tends to be piece-wise one-dimensional, whereas the configuration with no overlap (Config_3c) results in more complex two-dimensional fluid flow in the bedding plane. This need for increasing the number of nodes when decreasing the overlap between fractures is observed for both methods Inter1 and Inter2.Concerning the random fracture networks presented in Figures7 and 8(c)–8(f) show that both methods Inter1 and Inter2 can handle sparse configurations since the error in exiting flow rate Ex is smaller than 2% when (i) l is larger than or equal to 100 for both methods Inter1 and Inter2 when studying Config_4a (black lines with circles in Figures 8(c) and 8(d)) and (ii) l≥200 and l≥500 for method Inter1 and Inter2, respectively, when studying Config_5a (black lines with circles in Figures 8(e) and 8(f), respectively). However, increasing the percolation parameter p from 0.2 (Config_4a) to 0.25 (Config_4b) and 0.3 (Config_4c) for the sets of orthogonal fractures (Config_4) and from 0.5 (Config_5a) to 0.55 (Config_5b) for the sets of fractures with random orientations (Config_5) results in increasing the error Ex when using method Inter1 with Ex≈10% in Figure 8(c) (blue line with crosses and green line with stars) and 100% in Figure 8(e) (blue line with crosses). As illustrated in Figure 9, these large values of Ex are related to instabilities of the hydraulic head that are observed in the presence of small subfractures. In these cases, it is required to define the subfractures of the systems with the method Inter2 in order to provide a good representation of the fluid flow. Figures 8(d) and 8(f) show that using our modeling approach with the method Inter2 leads to values of Ex that are below or around 2% for all the random fracture networks presented in Figure 7.Figure 9
Hydraulic head distribution on the discretization nodes for the configurationConfig_4b (Figure 7(b)) using the methods (a) Inter1 and (b) Inter2 for defining the subfractures of the system. In (a), the illustration enlargement shows instabilities of the hydraulic head that are related to the presence of small subfractures, whereas these instabilities are not observed in (b).
(a)
Inter1
(b)
Inter2


## 5. Conclusions

In this work, we present the first attempt to develop three-dimensional hybrid analytical and numerical approaches for modeling fluid flow in fractured rocks. Although the results obtained in this study are satisfactory when defining the subfractures of the systems as pieces of the fracture borders (i.e., intersection methodInter2), additional developments and mathematical analysis are required for improving the stability of the developed solution. To this end, the next step of this work will focus on considering regularization methods that have been developed for different research fields using Fourier series as well [38, 39], which could help to stabilize the analytical solutions that are used at the subfracture scale in our models. This step is required before pushing further our efforts in terms of development of the method and comparison with existing modeling approaches.Future extensions will also be considered in order to increase the complexity of the fracture networks that can be handled with the method. Relying on the use of analytical solutions at the fracture scale with no source term (f=0 in Appendix A), the method is restricted to simplified fracture networks for which (i) the fractures are represented as rectangles and (ii) the intersection lines between the fractures are parallel to the fracture borders. The latter limitation comes from the technique that is used to couple the fractures, which relies on discretization nodes that must be located on the fracture borders. By using analytical solutions with f≠0 where this source term represents the flow entering or leaving the fracture through the intersection with another fracture, the restriction on the fracture intersection lines could be removed. Evaluating the feasibility of such an extension requires to (i) formulate the mathematical problem at the fracture scale with unknown discretized values of the hydraulic head gip located on both the fracture borders and fracture intersection lines, (ii) derive the corresponding outward fluid flow rates in order to define the new linear system Ax=b at the fracture-network scale, and (iii) implement these changes, including the definition of new discretization nodes, in the numerical method.Yet, we believe that the results presented in this work show the potential of three-dimensional hybrid analytical and numerical approaches, which are promising methods for providing efficient strategies to model flow processes in 3D fractured rocks. The strength of these methods relies on the discretization of the subfracture borders while no meshing is required in the fractures. This is a key advantage for future extensions to fluid and electric current flow applications that require taking into account the permeable surrounding rock matrix. Although these methods rely on simplified representations of the fractures and fracture networks, the need for efficient modeling approaches for data inversion and stochastic analyses forces us to explore various modeling options at the expense of the model complexity.



---

*Source: 3583817-2020-03-26.xml*