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

,


Introduction
Modeling fluid flow in the subsurface is required for numerous research fields and applications (e.g., [1][2][3][4][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 threedimensional domains (e.g., [7][8][9]).These representations are well suited to capturing the impact of fracture-network properties on large-scale observations (e.g., [10][11][12][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 fracturenetwork 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][19][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][22][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][25][26][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.

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 aperture b.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]): In equations ( 1), the hydraulic fracture transmissivity is given by T = ρgb 3 /ð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. 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.

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 Ω = fðx f , y f Þ: 0 ≤ x f ≤ L, 0 ≤ y f ≤ Hg where x f and y f are the local Cartesian coordinates defined along the subfracture borders, the corresponding boundary-value problem (BVP) is expressed as where Γ i ði = 1, 2, 3, 4Þ are the borders of the considered subfracture as shown in Figure 2.
As described in Appendix A, the solution hðx, yÞ of BVP (2) is expressed as with In expression (3), the functions g i describe the unknown hydraulic head h along the borders Γ i , which are discretized as shown in Figure 2. Discretizing each border Γ i into N i nodes results in discretizing the corresponding unknown function g i into N i discretized values denoted g p i ðp = 1, ⋯, N i Þ.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.

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: In expression (4), I m F represents the set of subfractures that share node m; J k ðx m k Þ is the average flow rate per unit length leaving fracture k at position x m k = ðx m k , y m k Þ, 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 = l m f /N m i with l m f the length of the fracture on which node m is located and N m i the number of discretized nodes on this fracture, which is determined from the following expression: N m i = min ð1, int ðl m f /Δ disc ÞÞ.In this expression, the function min ðx, yÞ returns the minimum value between x and y and subfractures with the method Inter2.The intersection line between the two main initial fractures is shown in red. 3 Geofluids the function int ðxÞ 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 g p i 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 g p i .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.

Method Analysis and Applications
The modeling approach described in Section 2 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 twointersecting 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 Figure 3 that correspond to a single fracture (Config_1a) and two intersecting fractures (Config_1b) in 10 × 10 × 10 m 3 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 hðx, 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.
Using the reference solutions described above, we define the relative error in the hydraulic head E h1 ðx, yÞ = |hðx, yÞ − h ref ðx, yÞ|/h ref ðx, yÞ × 100 with hðx, yÞ and h ref ðx, yÞ as the hydraulic heads at position ðx, yÞ computed with our hybrid approach and the reference solution, respectively.We also define E h2 ðxÞ as the average error at each position x of the relative errors E h1 computed along position y. Figure 4 shows these errors along the relative distance x L = 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 E h1 and E h2 with values below 2% and 1% when Δ disc ≤ 1 m and 0.5 m, respectively, except for some instabilities of E h1 that are observed when x L is located at the fracture extremities (i.e., x L close to 0 and 1 in Figure 4(a)).Conversely, when Δ disc = 5 m, large values of E h1 and E h2 are observed because only two discretization nodes are used to describe the hydraulic head distribution along the fracture borders.
The presence of instabilities of the hydraulic head shown at the fracture extremities in Figure 4(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 E x 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   4 Geofluids 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 E x for the configuration Config_1a.The results displayed in Figure 5(a) show that E x 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 E x , implying that we will use this value for the rest of this work.
In 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 E x for the configuration Config_1b (Figure 3(b)) considering the following definitions of ε: , where Δ m av f ract and Δ m av border 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 ε = Δ m av border 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 Section 2 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 m 3 .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)).

Geofluids
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.
In relation with previous studies describing fractured aquifers and reservoirs, we also consider the sets of random fracture networks Config_4 and Config_5 shown in Figure 7 with the percolation parameter p = ∑ N f i ð0:5l i Þ 3 /L 3 , with N f being the number of fractures, l i the length of fracture i, and L the domain size [24,30,36,37].In these 10 × 10 × 10 m 3 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 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.Our modeling approach is tested with these fracture networks by defining the relative errors E z and E x 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.

Results and Discussion
Figure 8 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, Figure 8(a) shows that the error in flow rate E z quickly decreases for the configurations Config_2 when increasing the characteristic discretization length l and reaches values below 2% when l ≥ l 1 and around 1% when l > l 1 with l 1 = 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 ≤ l 1 ), only one node is used on these smallest borders, and when Δ disc < 0:05 m (i.e., l > l 1 ), 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   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 E z is larger than 2% when l = l 1 .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 E z 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 Error in the exiting flow rate,  .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.9 Geofluids 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 Figures 7 and 8(c)-8(f) show that both methods Inter1 and Inter2 can handle sparse configurations since the error in exiting flow rate E x 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 E x when using method Inter1 with E x ≈ 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 E x 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 E x that are below or around 2% for all the random fracture networks presented in Figure 7.

Conclusions
In this work, we present the first attempt to develop threedimensional 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 method Inter2), 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  10 Geofluids the mathematical problem at the fracture scale with unknown discretized values of the hydraulic head g p i 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.

Appendix A. Analytical Expression for the Hydraulic Head
The solution hðx f , y f Þ for the boundary-value problem (BVP) with 0 ≤ x f ≤ L and 0 ≤ y f ≤ H, can be expressed as [40,41] with Considering f ðx f , y f Þ = 0, and using the properties of the Dirac delta function δ provided in Butkovskiy [40] and of trigonometric and hyperbolic functions provided in Gradshtein and Zwillinger [42], hðx, yÞ is expressed as x f , i = 3, 4: ( ðA:2bÞ

B. Analytical Expressions for the Outward Fluid Flow Rates
Consider a subfracture of length L, width H, and constant transmissivity T, the outward averaged flow rates per unit length J i on subfracture border Γ i ði = 1, 2, 3, 4Þ are defined as ðA:1cÞ 11 Geofluids where, as described in Section 3, ε is the distance parameter that is used to reduce the instabilities observed at the extremities of the Fourier series.
Using expression (3) and discretizing the integrals related to functions g i ði = 1, 2, 3, 4Þ with N x and N y nodes in x f -and y f -directions, respectively, the derivatives of hydraulic head required in (B.1) are expressed as with g i ðη p Þ = g p i ði = 1, ⋯, 4Þ, Δη = H/N y , and Δξ = L/N x .

Figure 1 :
Figure 1: Example of (a) two intersecting fractures that have been divided into (b) 8 subfractures with the method Inter1 and (c) 4 subfractures with the method Inter2.The intersection line between the two main initial fractures is shown in red.

Figure 2 :
Figure 2: Illustration of a subfracture having length L, width H, local coordinates x f and y f , and discretized borders Γ i ði = 1, 2, 3, 4Þ where the discretization nodes are represented with crosses.

𝛥 disc = 5 𝛥 disc = 1 𝛥 disc = 0. 5 𝛥Figure 4 :
Figure 4: Errors in hydraulic head along the relative fracture length x L = 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 E h1 at each position ðx, yÞ and (b) the relative error E h2 averaged over position y.The dashed and dotted horizontal black lines represent errors of 1 and 2%, respectively.

Figure 5 :
Figure 5: Error in the exiting flow rate E x (%) 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.
Figure 8 shows E z for the sets of fracture networks Config_2 and Config_3 and E x 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 6 :
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.

7 models
Geofluids Config_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 E z is represented in Figure 8(b) with full and dashed lines, respectively.

Figure 7 :
Figure 7: Random fracture networks considering (a-c) orthogonal families of fractures and (d, e) random orientations for various values of the percolation parameter p.

Figure 9 :
Figure 9: Hydraulic head distribution on the discretization nodes for the configuration Config_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).