Impact of Fracture Topology on the Fluid Flow Behavior of Naturally Fractured Reservoirs

Abstract

Fluid flow modeling of naturally fractured reservoirs remains a challenge because of the complex nature of fracture systems controlled by various chemical and physical phenomena. A discrete fracture network (DFN) model represents an approach to capturing the relationship of fractures in a fracture system. Topology represents the connectivity aspect of the fracture planes, which have a fundamental role in flow simulation in geomaterials involving fractures and the rock matrix. Therefore, one of the most-used methods to treat fractured reservoirs is the double porosity-double permeability model. This approach requires the shape factor calculation, a key parameter used to determine the effects of coupled fracture-matrix fluid flow on the mass transfer between different domains. This paper presents a numerical investigation that aimed to evaluate the impact of fracture topology on the shape factor and equivalent permeability through hydraulic connectivity (f). This study was based on numerical simulations of flow performed in discrete fracture network (DFN) models embedded in finite element meshes (FEM). Modeled cases represent four hypothetical examples of fractured media and three real scenarios extracted from a Brazilian pre-salt carbonate reservoir model. We have compared the results of the numerical simulations with data obtained using Oda’s analytical model and Oda’s correction approach, considering the hydraulic connectivity f. The simulations showed that the equivalent permeability and the shape factor are strongly influenced by the hydraulic connectivity (f) in synthetic scenarios for X and Y-node topological patterns, which showed the higher value for f (0.81) and more expressive values for upscaled permeability (kx-node = 0.1151 and ky-node = 0.1153) and shape factor (25.6 and 14.5), respectively. We have shown that the analytical methods are not efficient for estimating the equivalent permeability of the fractured medium, including when these methods were corrected using topological aspects.

1. Introduction

Naturally fractured reservoirs (NFRs) are composed of various lithologies such as shales, sandstones, carbonates, and igneous rocks and correspond to over 30% of the global production [1]. Within this group, fractured carbonate reservoirs represent an important part of the world’s oil and gas reserves (e.g., the Campos and Santos Basins in Brazil and the Kwanza Basin in Angola) [2,3]. Characterization and numerical simulation of naturally fractured reservoirs represent a challenge because of their complex evolution (diagenesis, geomechanics) and the effect of the coupled fracture-matrix relationship on the fluid flow [4,5,6,7,8].

Mathematical modeling of fluid flow and transport in a fractured reservoir also represents a challenging task. First, the acquisition of fracture systems attributes (e.g., aperture, length, orientation, spatial distribution, and topology) is complex because of the need to integrate various scales and limited sampling through subsurface datasets (e.g., seismic, well logs, and cores) [9,10,11]. The scale of the fracture systems is an important aspect to be considered in the fluid flow analysis of naturally fractured reservoirs because it affects the connectivity of the fracture network (spacing, orientation, clustering, and fracture intensity), as discussed by Nelson [12], Laubach [13], Berkowitz et al. [14], Roy et al. [15], Laubach et al. [16], Sahu and Roy [17], and Silva et al. [18], among others. Second, the main challenges are posed by the anisotropic nature of fracture systems controlled by stress distribution and the rock heterogeneities [19] and the dual porosity phenomenon when the rock matrix and the surrounding fractures behave as separated media with variable grades of interaction [20]. Large fractures represent a crucial aspect in fractured systems because it is often necessary to model them explicitly due to their singular effect in numerical models [21].

There are different ways of modeling an NFR. One is the approach of equivalent continuous medium, which assumes that the fluid flow in the geological medium, including fractures, can be simulated in a representation of a porous medium with the classical equations of flow and transport Another approach is the discrete fracture network model (DFN), which considers the representation of fractures in the model, and the flow simulation considers their influences. The main differences to flow simulation between dual-porosity and discrete fracture models are related to time flow behavior during the early phases of simulation because the saturation is modeled with different efficiency [22]. A third approach to simulate an NFR is called hybrid-dimensional modeling and involves the discretization of major faults and fractures, with the matrix and small fractures represented as a continuous medium [23]. The main objective is to facilitate mesh generation. DFN models allow better prediction of fluid flow and transport phenomena than continuum models [24,25]. An advantage of the discrete fracture technique is that it can explicitly consider the effects of discrete fractures on reservoir fluid flow [26,27,28]. Thus, DFN models have become popular due their practical applications [28]. The fluid flow through NFRs is highly heterogeneous, as flow develops through the matrix and fractures because to the hydraulic head gradient [29,30,31]. Fluid flow in the rock matrix is limited by the matrix permeability and if the fracture network allows, pressure gradients will form in each block [32,33]. Thus, the permeability of the NFRs is not an intrinsic property but changes strongly with the sampling and scale of characterization [34]. Embedding fractures as discontinuities in finite element mesh (FEM) can adequately represent the behavior of natural fractures in a hydraulic [18] and hydromechanical numerical analysis [35,36]. This technique allows us to estimate equivalent permeability in naturally fractured media. The advantage of using DFN models to compute the equivalent permeability of fractured rocks through numerical analysis is linked to the upscaling approach, which allows reaching computational efficiency and generating reliable input data for double porosity-permeability models. The upscaling process is achieved with the replacement of the fracture network with homogeneous grid blocks, and the average flow can be adequately estimated [18,36,37,38,39,40,41].

The main objective of this work is to investigate the impact of fracture systems geometry, based on analysis of topology, in the fluid flow in naturally fractured media. Two fracture systems can exhibit the same geometrical elements regarding orientation and trace length but present different node (where fracture planes terminates or cross other fractures) arrangements [42]. The fracture topology can have different impacts on the flow [41,42,43,44]. The hydraulic connectivity (f) is related to the fracture topology. This variable allows obtaining an analytical value of equivalent permeability for the fractured medium using Oda’s [45] methodology, which considers discontinuous rock masses as homogeneous, anisotropic porous media.

The double porosity approach is a computationally efficient and usually adopted method to simulate the flow in the fracture-matrix system. The formulation of double porosity-permeability involves the determination of fluid transfer between the matrix and the fracture, related to the shape factor as previously mentioned, and this approach has been evaluated by several authors [46,47,48,49,50,51,52,53]. Although several strategies have been evaluated to enhance the predictions of fracture-matrix exchange in double porosity-permeability models, most commercial formulations are based on the initial model of fractured reservoirs as proposed by Warren and Root [46,54]. For irregular fracture networks, which are more realistic, there is not a unique geometric method to estimate the shape factor. Performing numerical simulation of models with discrete fractures aims to predict the flow behavior [55]. Thus, we have performed numerical simulations to compare the effect in the flow patterns regarding the application of the shape factor. We have used hypothetical and real scenarios with different fracture topologies (I-nodes, Y-nodes, and X-nodes) to provide a realistic investigation of the problem.

Our results showed the impact created by the arrangements of fracture nodes on the fluid transfer term between the rock matrix and the fracture network, which is directly related to the capacity of the fluid flow in the system, the shape factor, and the equivalent permeability. The contribution of this paper is related to the proposition of an integrated methodology that considers the adoption of topological pattern analysis for computing the hydraulic connectivity number, regarding the influence of this parameter for the numerical and analytical upscaling process of equivalent permeability and dual porosity dual permeability shape factor for naturally fractured reservoirs.

2. Mathematical Formulation

The methodology adopted used DFN-based models to obtain equivalent permeability fields and the shape factor and compared with solutions obtained for flow properties using Oda’s model. This procedure considered the hydraulic connectivity (f) to estimate the influence of each topological configuration. The shape factor was also analyzed considering the hydraulic connectivity f from the topological analysis. The mathematical form of the flow equation for double porosity-permeability systems is presented in Equation (1), as proposed by Lewis and Schrefler [56].

As mentioned above, the shape factor and equivalent permeability are input values for numerical simulation of the double porosity-permeability model. In Equation (1), the term (α) refers to the fluid transfer term between the fractures and the rock matrix in a naturally fractured medium and is related to the shape factor. Permeability is defined within the source term in Equation (1) and is presented in more detail in Equations (2) and (3) (k_i).

(1)${\mathsf{\varphi }}_{\mathrm{i}}\frac{\partial {\mathsf{\rho }}_{\mathrm{f} \mathrm{i}}}{\partial \mathrm{t}}+\nabla ·\left[{\mathsf{\rho }}_{\mathrm{f} \mathrm{i}}{q}_i\right]+{\left(-1\right)}^{\mathrm{i}+1}\mathsf{\alpha }=0\hspace{1em}\hspace{1em}\mathrm{i}=1,2,$

(2)$q_i=-K_i\nabla {\mathrm{P}}_{\mathrm{i}},$

(3)$K_i=\frac{k_i}{\mathsf{\mu }},$

where ϕ is the porosity, ρ_f is the fluid density, q is the fluid flux vector, P is the pressure and μ is the fluid viscosity.

Below, we describe each of the stages of the methodology.

2.1. Topological Characterization of Fracture Network

An NFR may involve several fracture sets, which may or may not intersect each other, and fractures can be described by a variety of geometrical attributes [28]. The geometry of the fractures and the relationships between individual fractures or fracture sets control the flow properties [41]. Regarding the importance of fracture connections, the term topology refers to the spatial relationships between fracture planes, and the topological analysis can improve the analytical estimation of effective permeability in fracture networks [41]. Topology defines the connections between fractures in a discretized fracture network model as node counting and branch analysis. The proportions of abutting Y-, crossing X-, and isolated I-nodes provide topological classification and node type counting (Figure 1). Topology analysis can determine the dimensionless intensity, frequency, orientation, and length of fractures. Sampling proportions of various types of nodes provides a topological classification that can estimate the conductivity effect of the network [42] (Figure 1).

Sævik and Nixon [41] presented an innovative analytical method for the estimation of effective permeability using topological parameters. They focused on analytical upscaling to create a functional relationship between the effective permeability and the fracture parameters.

The hydraulic connectivity f is the ratio between the actual effective permeability of the network (upscaling) and the theoretical upper bound as defined by Oda [45]. In this framework, the connectivity f considers the value of 1 (f = 1) in an idealized case of fractures defined by two parallel planes, a classical representation of rock fractures in numerical models [28]. For other idealized fracture geometries, f is a function of the fracture intensity, defined as the total trace length per sampling area [41].

The value of hydraulic connectivity f is computed by conducting a numerical simulation of fluid injection through the grid. The equivalent permeability obtained by numerical simulation was also processed through Oda’s analytical method to define f. The topological connectivity parameters were calculated and plotted with parameter f. We have assumed that the hydraulic connectivity f is zero below the percolation threshold.

Using the maximum operator, the equation can be written as follows:

(4)$f=max(0,\frac{2.94·\left(4n_x+2n_y\right)}{4n_x+2n_y+n_i}-2.13),$

where n_x is the number of nodes X, n_y is the number of nodes Y and n_i is the number of I type nodes.

2.2. Discrete Fracture Model and Finite Element with Embedded Discontinuities

We used the mathematical development of the continuous approximation of embedded discontinuities proposed by Beserra et al. [35], which was also applied by Silva et al. [18]. This approach was implemented in the CODEBRIGHT in-house software [57,58], a finite element-based code to perform numerical analysis of fluid flow in deformable porous media in a fully coupled scheme. The embedding process, which considers the triangular element with three nodes from the domain (length, band, and width), is divided into two parts, isolating node 1 from nodes 2 and 3, as shown in Figure 2a. Darcy’s flow into the discontinuity is decomposed in a tangential and normal component as is shown in Figure 2b.

The formulation adopted in the present work considers that the fluid flow in a fracture occurs in its direction. Therefore, we considered only the tangential component. Darcy’s law can be written for continuum (q_Ω) as follows:

(5)$q_{\mathsf{\Omega }}=-K_{\mathsf{\Omega }}\nabla \mathrm{p},$

where K_Ω is the permeability tensor of the continuous part, and Darcy’s law for fracture (q_s) is

(6)$q_{\mathrm{s}}=-{\mathrm{K}}_{\mathrm{s}}\nabla \mathrm{p},$

where K_s = k_s/μ.

The finite elements approach considers that the flow at the fracture will occur over the entire element thickness. As presented by Silva et al. [18], to ensure that the element transmissivity in the fracture direction is the same as the transmissivity of the embedded fracture, the fracture permeability K_s is multiplied by the ratio h/l_e. The resulting Equation is:

(7)$q_s=-\frac{h}{l_e}{\mathrm{K}}_{\mathrm{s}}t\otimes t\nabla \mathrm{p},$

where t is the discontinuity direction, h is the thickness of incorporated discontinuity, and l_e is the element size.

Because a fracture in a porous medium represents a preferred path for flow toward discontinuity, it induces anisotropy in the direction of the fracture throughout the medium. As a result, an effective permeability tensor for the element is given by:

(8)$q=-K_{\mathrm{ef}}\nabla \mathrm{p}$

(9)$K_{\mathrm{ef}}=\left(K_{\mathsf{\Omega }}+\frac{h}{l_e}{K}_{s}t\otimes t\right)$

Two different permeability variation laws are adopted, one for the continuum portion and the other for the fracture. Estimation of intrinsic fracture permeability k_s can be performed through models that consider fractures with rough walls or filled with proppant material. We used the law of parallel plates in the models built for this research [35,59].

(10)$k_s=k_o{\left(\frac{h}{h_o}\right)}^2$

where k_o is the permeability in the continuum portion, h is the fracture thickness, and h_o the initial fracture aperture.

For the rock matrix, we adopted the isotropic permeability law, described as follows:

(11)$k=k_o{\left(\frac{\varphi }{{\varphi }_o}\right)}^A$

where ϕ_o is the initial porosity, ϕ is the porosity calculated after perturbations induced by rock deformation, and A is a parameter to be determined experimentally [36,60].

The continuity equation that governs the fluid flow in the porous media, considering a rigid medium and without gravitational effect, is given by:

(12)$\mathsf{\varphi }\frac{\partial {\mathsf{\rho }}_{\mathrm{f}}}{\partial \mathrm{t}}+\nabla ·\left[{\mathsf{\rho }}_{\mathrm{f}}q\right]=0$

where ϕ is the host rock porosity and q is the fluid flux vector given by Equation (8).

The fluid density is related to pore pressure according to the following expression:

(13)${\mathsf{\rho }}_{\mathrm{f}}={\mathsf{\rho }}_{\mathrm{f}0}{\mathrm{e}}^{\mathsf{\beta }\left(\mathrm{p}-{\mathrm{p}}_0\right)}$

where ρ_f0 is the reference density, β is the fluid compressibility, p is the pressure, and p₀ is the reference pore pressure.

2.3. Calculation of Shape Factor for Dual-Porosity—Permeability Models

The finite element method with embedded discontinuities allowed monitoring the fracture flow rate by applying the fluid flow Equation in the rock matrix to derive the shape factor that best corresponds to the simulation results.

Assuming that the matrix is a homogeneous and isotropic medium, ϕ_m represents the porosity, β_m the compressibility, k_m the permeability tensor, and μ the viscosity, a general expression of the Equation for matrix pressure, p_m(x,t), is given by:

(14)$\{\begin{array}{c}{\varphi }_m{\beta }_m\frac{\partial p_m\left(x,t\right)}{\partial t}+\nabla ·q_{\mathrm{m}}=0 em \Omega \times \left(0,T\right)\\ q_{\mathrm{m}}=-\frac{k_m}{\mu }\nabla p_m\left(\mathit{x},t\right)\\ p_m\left(x,t\right)=p_f \text{ for } {\Gamma }_{mf}\times \left(0,T\right)\\ q_{\mathrm{m}}·n=0 \text{ for }{\Gamma }_N\times \left(0,T\right)\\ p_m\left(x,\mathrm{t}\right)=p_{m,0}\text{ in }\Omega \times \left(t=0\right)\end{array}$

Considering the domain border Ω given by $\partial \Omega ={\Gamma }_{mf}\cup {\Gamma }_N$ where Γ_mf represents the matrix-fracture interface, Γ_N is the boundary on which the null flow is imposed, and n is the normal vector of Γ_N outside the domain Ω, p_f is the fracture pressure, and p_m,0 is the initial matrix pressure.

Zimmerman et al. [59] developed two equivalent expressions to calculate the shape factor, one based on the mean matrix pressure and the other with flow defined by $\partial \Omega$. Thus, the following expression is adopted to introduce the shape factor α(t):

(15)$\nabla p_m\left(x,t\right)·n\approx \frac{\overline{p_m}\left(t\right)-p_f}{\delta t}$

where $\overline{p_m}\left(t\right)$ is the measure of p_m(x,t) on domain Ω:

(16)$\overline{p_m}\left(t\right)=\frac{1}{\left|\Omega \right|}{{\displaystyle \int }}_{\mathsf{\Omega }}{p}_m\left(x,\mathrm{t}\right)d\Omega$

The parameter $\delta t$ in Equation (15) measures the penetration depth of perturbation into p_m(x,t) caused by p_f. Thus, the expression to calculate the shape factor is given by:

(17)$\alpha \left(t\right)=\frac{\left|{\Gamma }_{mf}\right|}{\left|\Omega \right|\delta t}$

where |Γ_mf| denotes the area Γ_mf, and α(t) is the transient shape factor.

We have adopted an expression to obtain the shape factor as a function of mean pressure $\overline{p_m}$, using the Gauss theorem and the balanced Equation (4), described as follows:

(18)$\alpha \left(t\right)=\frac{{\varphi }_m{\beta }_m\mu }{k_m\left(\overline{p_m}\left(t\right)-p_f\right)}\frac{\partial \overline{p_m}\left(t\right)}{\partial t}$

In the limit t → ∞, the stationary value of such a parameter is obtained.

(19)${\alpha }_{\infty }=\underset{t\to \infty }{lim} \alpha \left(t\right)$

As α_∞ is a geometric parameter, it does not depend on p_f or the hydrodynamic properties of the matrix.

2.4. Oda Tensor Approach for Equivalent Permeability Estimation

Oda [45] proposed an analytical method to compute equivalent fracture permeability, which is widely used in DFN upscaling due to its efficiency [61,62]. For a grid cell with known fracture areas, A_k, and transmissivities, T_k, obtained from the DFN model, an empirical tensor that considers only fracture geometry can be calculated as follows:

(20)$F=\frac{1}{V}{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{N}}{f}_k{A}_k{T}_{k}n\otimes n$

where F is the fracture tensor, V is the grid cell volume, N is the total number of fractures in a grid cell, f_k is the hydraulic connectivity for fracture k (generally assumed to equal to 1), A_k is the area of fracture k, T_k is the transmissivity of fracture k, and n is the components of a unitary normal to the fracture.

Oda’s permeability tensor is derived from F by assuming that F expresses fracture flow as a vector along the fracture’s unit normal. Assuming that fractures are impermeable in a parallel direction to their unit normal F must be rotated into planes of permeability [63].

(21)$k=\frac{1}{12}\left(\mathrm{tr} F I-F\right)$

where k is the permeability tensor, F is the fracture tensor, I is the identity matrix, and tr F is the sum of the main diagonal of the fracture tensor. According to Elfeel [64], a limitation to the application of Equation (21) is that fractures of any length will contribute to the permeability value, including if fractures do not allow percolation in a cell.

3. Modeling of Studied Cases

The characteristics of models, four synthetic cases, and three real scenarios of a Brazilian pre-salt reservoir, used to perform the proposed methodology are described below.

3.1. Artificial Fractured Media

The synthetic 2D sections analyzed in this work were created based on scenarios presented and tested by Silva et al. [18]. We have adopted changes in their models regarding the dimensions and fracture characteristics as the apertures. The scenarios consider four different topological patterns with different intensities and distribution of fractures that allowed the analyses of the impact of the topology on the shape factor and validating the results they present for the equivalent permeability (Figure 3).

The quantification of topological fracture node types for each scenario (Figure 3) is presented in Table 1.

For all cases, the fracture intensity was calculated according to the method proposed by Dershowitz and Einstein [24], which defines measures of fracture abundance by Pxy, where x is the dimension of the sampling region and y the dimension of the measured feature (e.g., Sanderson and Nixon [42]).

We used the P22 sampling technique (area of fractures per sampling area) to calculate the relationship between the topology of the fractures and the equivalent permeability. The dimensionless intensity calculated from the P22 analysis is given by $P_{22}=N_L{L}_C^2/A=P_{21}·L_c$, and it has dimensions [L⁻⁰] where L_c is the characteristic length and is defined as the arithmetic mean of the line lengths. We used a MATLAB^® script to build the FEM (parameterized meshes). After loading the fracture coordinates (endpoints) and apertures, the script allows the extraction of the fracture permeability using the parallel plate theory (Equation (10)). We consider the matrix as isotropic and of low permeability, which simulates a tight reservoir scenario.

The simulation of hydraulic phenomena consisted of the injection of a single fluid-phase (water) on the left side of the models and fluid production on the right side of the domain. The fluid injection pressure was 55 Mpa, and the production pressure is 54.9 Mpa, leading to a pressure gradient of 0.1 MPa. The simulation was also conducted for all scenarios considering a vertical flow, using the same pressure values (Figure 4).

For the simulations, we considered a single-phase flow governed by a nonlinear form of Darcy’s law, the equivalent permeabilities computed from DFNs, and a flow-based upscaling approach.

Thus, the numerical simulation was conducted to obtain the shape factor. In the FEM used for simulations, the fractures were considered as flow lines. The shape factor was calculated numerically using the discrete fracture model for all flow states until reaching its pseudo-stationary value.

3.2. Real Fractured Reservoir Case: Brazilian Pre-Salt

The real case analyzed in our study is adapted from Falcão et al. [36]. They developed several DFN models for a carbonate reservoir from the Brazilian pre-salt interval (Figure 5). Figure 5 shows the model used to extract fracture distribution examples (bidimensional sections) used in this work. The Santos Basin pre-salt cluster is located in ultra-deep waters, 2000 m approximately, and 290 km offshore from the Rio de Janeiro coast. This cluster represents a prominent regional sequence covered by a continuous evaporitic sequence. The reservoirs are mainly composed of carbonate rocks, laminites, stromatolites, coquinas, and spherulites. Pre-salt development presents many challenges like water depths, reservoir heterogeneities, oil with the content of CO₂ in the gas phase, and drilling challenges because of variable salt layer [65].

Figure 6 shows the three bidimensional scenarios, horizontal sections (slices), extracted from the tridimensional cell from the pre-salt reservoir model shown in Figure 5c. Each 2D section (200 × 200 m) presents a different fracture intensity and spatial arrangement. We have adopted the fracture parameters and coordinates used to compose the 2D sections after the original work of Falcão et al. [36]. The first 2D model, section 1, possesses 54 fracture planes (Figure 6a), the second, section 2, possesses 22 fractures (Figure 6b), and the third case, section 3, 75 fractures (Figure 6c). The building of the FEM created from the sections shown in Figure 6 followed the same procedures used to build the FEM from the artificial models shown in Figure 3 and Figure 5.

Table 2 shows topological pattern values of each section shown in Figure 6. These values are used in Equation (4) for the calculation of the hydraulic connectivity value f.

4. Results

4.1. Artificial Scenarios

Figure 7 shows the results for the three artificial cases regarding the equivalent permeability calculated with the horizontal flow simulation (direction x). Oda’s method overestimates the equivalent permeability because it assumes that all the fracture planes cross the entire sample. Permeability calculated using the Oda method with the hydraulic connectivity f reached a better result, similar to the numerical simulation for the sections with Y- and X-nodes. The calculation of permeability with (f) for the I-nodes case resulted in negative values (not shown in the graph) because the values of n_x and n_y are zero, which resulted in inconsistent outputs, and the 0 value is assumed.

It indicates that Oda’s method cannot compute the equivalent permeability for the situation of topology dominated by I-nodes, which has no connection between fracture planes. The corrected Oda’s method presented reliable results for the cases dominated by X- and Y-nodes. The same aspect was demonstrated by Madden [66], that the geometric mean of heterogeneous bodies could be used to predict the physical properties of a homogenized matrix.

Figure 8 shows the results obtained for the equivalent permeability values obtained through the vertical flow simulation (y-direction). Results are close to those obtained for the horizontal flow (Figure 7) and agree with Silva et al. (2021). This similarity occurs because the arrangement does not affect the flow pattern in the two situations (Figure 7 and Figure 8).

Figure 9 presents the shape factor behavior during the simulation for the synthetic scenarios. The significant impact of the I-node-dominant topology is shown, and the related curve reaches the pseudo-stationary state first, followed by the X-node and Y-node curves (Figure 9). This effect demonstrates that the lack of fracture connections dramatically affects the fluid transmissibility through the matrix.

The shape factor is defined when the curves flatten parallel to the X-axis. The building of synthetic cases has considered similar dimensionless fracture intensities based on the P22 distribution. This condition is important to assure that similar fracture areas with different element distributions would be compared.

Figure 10 presents the dimensionless fracture intensity P22 percentage defined for each synthetic case (Figure 10a) and the final shape factor value calculated for each case (after reaching the pseudo-stationary state) (Figure 10b). The case dominated by the X-node pattern shows the highest value of shape factor as expected because of the higher connectivity relationships. A similar response is observed for Y-node dominated topological patterns, while for combined (XYI) and dominated I-nodes scenarios, the values are low. It leads to the evidence that the fracture crossing predominance induces a more considerable matrix-fracture transfer and can impact the fluid flow for the dual-porosity dual permeability media.

Figure 11 shows the result of a cross-correlation between the computed equivalent permeability and shape factor for each case. The patterns dominated by X-nodes and Y-nodes show slightly higher equivalent permeability, which demonstrates the negative impact caused by I-nodes in the X-Y-I pattern (Figure 11), and that the shape factor, regardless of the fracture intensity, affects the permeability anisotropy [67].

Table 3 shows the properties of the synthetic cases and results obtained through the simulation: the equivalent permeability, factor f, and the shape factor. It demonstrated that the lower shape factor values resulted from the influence of the I-node terminations.

4.2. Numerical Simulation of Real Fractured Reservoir Scenarios of Brazilian Pre-Salt Interval

Figure 12 and Figure 13 show the results of equivalent permeability computed through Oda’s method, and Oda’s method corrected with hydraulic connectivity and numerical simulation for horizontal and vertical flow for the three real cases (Figure 6). For the horizontal flow situation, Oda’s method yielded higher values of equivalent permeability, and it overestimated the numerical simulation results, which showed slightly lower values (Figure 12). Application of the Oda method corrected with the hydraulic connectivity f resulted in the lowest values for the equivalent permeability, compared to the numerical simulation. This effect results from the influence of I-node fractures, and the opposite is observed for Oda’s analytical solution (Figure 12). De Dreuzy [68] demonstrated a significant coupling between flow heterogeneities at the fracture scale and flow heterogeneities at the network scale that cannot be solved through the hydraulic connectivity approach, as observed in real cases presented here. For section 2, the value of hydraulic connectivity f was 0 due to the high number of I-node fractures compared to the X- and Y-node types (Figure 6b). Thus, Equation (4) becomes negative, which implies a lack of coherency.

Figure 13 presents the results for the calculation of permeability regarding the vertical flow scenario. The results showed higher values, in certain cases, more than the double extracted with each method, compared to the horizontal flow scenario (Figure 12). It occurs because the fracture planes are aligned with the flow direction. Values computed for sections 1 and 3 show a significant difference between the numerical and analytical results obtained with Oda’s model. As observed with the synthetic cases, this effect occurs because this method assumes that the fracture planes cross the entire model area. After the application of the correction factor to the topological aspect, the result becomes close to the numerical solution. Again, this result shows the influence of the I-node fractures in the transmissibility phenomena.

Figure 14 and Figure 15 show the cross-correlation of the computed equivalent permeability and the hydraulic connectivity f calculated with the topological analysis for the horizontal and vertical flow scenarios, respectively, for the three real cases treated. For the horizontal flow scenario, results obtained with Oda’s method are the closest to the result reached with numerical simulation. For this specific scenario, the hydraulic connectivity does not help to achieve a reliable solution. However, the value of hydraulic connectivity f is efficient to demonstrate the impact of the topological I-nodes in the three real cases analyzed.

For the vertical flow scenario, sections 1 and 3 show a similar correlation for the factors analyzed, and section 2 presents a remarkable difference relate to the connectivity effect (Figure 15). As observed for the horizontal flow scenario, the hydraulic connectivity f has a similar pattern because it does not depend on the direction of the flow in the fractured medium.

Figure 16 shows the shape factor calculated as a function of time for the three real 2D sections. The curves show the influence of fracture density on fluid transfer between the fracture network and the matrix. For section 2, which has only 23 fractures, it takes longer to reach a pseudo-stationary state. As observed above, the lack of coherence for the connectivity factor for section 2 prevents a better comparison of correlation with other parameters.

Figure 17 presents a cross-correlation, a comparison between the dimensionless fracture intensity P22 (the relationship between the area of fractures over the area of the section) Figure 17a, and the shape factor related to the topology patterns of fracture networks for the three 2D sections derived from the pre-salt reservoir model Figure 17b. The section 3 presents the higher fracture intensity and the lower shape factor (Figure 17). The section 2, dominated by I-node endings, shows the higher shape factor and the lower fracture intensity.

Figure 18 shows cross-correlation between the shape factor value at its stationary state and the calculated hydraulic connectivity f for each 2D section. Section 3 presents lower values for these two parameters, due to the influence of I-node fracture endings. However, section 3 presents the highest fracture intensity, as shown in Figure 17. Section 2 presents zero value for the shape factor due to the lack of connectivity caused by the dominance of I-node fractures, and the lack of coherence prevented correlation. Section 1 presents a higher connectivity shape value and, relatively, the higher shape factor of the three sections.

Figure 19 shows the correlation between computed equivalent permeability and the shape factor calculated from the topological scenarios with X- and Y-node types. Section 3 presents a lower shape factor but higher equivalent permeability. Section 2 shows markedly higher values of shape factor but slightly lower equivalent permeability values than section 3. Section 2 also presents the highest shape factor but relatively lower equivalent permeability to sections 1 and 3.

Table 4 shows the properties of the three 2D sections extracted from the pre-salt reservoir model and results obtained through the simulation: the equivalent permeability, factor f, and the shape factor. This integration highlights the relationship between the topological patterns and the resulting equivalent permeability. It also points to the importance of defining both the number of connections between fracture planes and the number of I-node types present in the DFN because it controls the permeability effect and the interaction of fluid between fractures and rock matrix.

We present a ternary diagram in Figure 20a that shows the percentage of each topology pattern for the three sections studied. Figure 20b shows a bar chart with the value for the hydraulic connectivity factor for the same sections. The hydraulic connectivity f tends to be higher in the sections with a higher percentage of X- and Y-nodes. However, it does not indicate higher equivalent permeability fields, and the definition of topological aspects may be essential to constraining hydraulic simulations in fractured media.

The results showed the main effects of topology variations in upscaled permeability tensor estimation and shape factor, considering a low and fixed matrix permeability (e.g., tight reservoir). This work provided an analysis of fracture patterns influence (type and degree of connections) in the cell scale used in reservoir simulation.

5. Discussion

There are many advancements since the early works on the numerical modeling of geological scenarios. Jing and Stephansson [44] showed that applying a mathematical formulation to represent fractured rocks cannot be arbitrary, and it must relate to the topological properties of the fracture network. Haddad et al. [69] showed that the heterogeneity of fractured reservoirs profoundly affects the rate of mass transfer between the fracture system and matrix blocks. One of the heterogeneity sources is related to the variable fracture intensity.

In this study, we successfully used embedded discontinuities in FEM-based models, regarding real and artificial scenarios, of fractured porous media to execute flow simulations. A similar approach was also successfully used for real cases in Brazilian pre-salt reservoirs by Falcão et al. [36].

Performing numerical simulation with the incorporated discontinuities has demonstrated reliability for fluid flow simulation in fractured media, including coupled simulations [70,71,72,73,74,75,76]. This method presents an alternative to the classical techniques based on the use of interface elements. It also represents an alternative for analysis that requires the repeated rewiring of the geometry, as observed in discrete fracture methods in which fractures are represented explicitly in the model [77,78]. The approach allows representing fractures with lines and surfaces [79]. This technique also presents the advantage of reducing the computational costs and simplifies the discretization process of fractured media [35] and directly includes the fracture contribution into the matrix without any spatial portioning.

The artificial and real case models treated here were built with the FEM method, and the flow simulations performed allow calculating the equivalent permeability and the shape factor for all the scenarios. The 2D real scenarios were built from a cell extracted from a larger DFN model of a reservoir from the Brazilian pre-salt interval [36]. Three horizontal slices were acquired from a 3D cell and provided three bidimensional sections with proper fracture sets. The simulations allowed us to evaluate the impact of fracture topology and highlight aspects of the computation of equivalent permeability and shape factor in the fractured medium. Fracture topology handles the generalization of spatial relationships as the connectivity and continuity of fracture planes. It refers to properties that are unchanged by the continuous transformation of the space in which the fractures are embedded, the deformation (strain) of the geometrical properties (lengths and angles), a concept familiar in many areas of geology [18,42]. Topological analysis of NFRs is fundamental to predict the hydraulic connectivity f, which involves the fracture network and the rock matrix in classical dual-porosity problems. Besides the numerical simulation, we used Oda’s analytical method, and the correction applied to this method according to Sævik and Nixon [41]. The correction considers the fracture system connectivity and fractures characteristics using the quantification of topological fracture node types. Although these authors do not compare results with Oda’s original method, several authors have made advances in this analytical method and shown the need to consider other representative parameters to represent the permeability tensor of fractured media more realistically [78,80,81]. Fumagalli et al. [78] implemented the embedded discrete fractures technique and compared the simulations with Oda’s method used in a commercial simulator. Their results showed the accuracy of the embedded fracture models compared to the analytical method.

Considering the correction to Oda’s method provided by the topological analysis, the hydraulic connectivity f is used to represent a new parameter in the general equation of Oda’s tensor. Results obtained through simulation for the synthetic cases with X- and Y-nodes were similar to the values obtained with the corrected Oda method. This outcome demonstrates the importance of analyzing the hydraulic connectivity f to estimate equivalent permeability [43,82].

Silva et al. [18] presented a numerical tool for fracture modeling, which was tested in synthetic bidimensional models of fractured media. They analyzed the capacity of FEM models with embedded discontinuities to calculate upscaled equivalent permeability and the influence of topological patterns in the flow. These authors analyzed the influence of topological patterns in the equivalent permeability. They concluded that the fracture intensity, connectivity, and topology pattern play an important role in the upscaled permeability. These authors demonstrated the reduced equivalent permeability of fracture networks dominated by I-nodes, due to the lack of connectivity, compared to networks with X- and Y-node patterns, as was also found by the present study.

Lahiri [83] performed the calculation of the effective permeability of a fracture network and used topological connectivity through analytical expressions, considering fractal and multifractal distributions. The author concluded that the permeability of fracture systems increases rapidly with the growth in length of the branch segments. Longer branch segments in fracture networks play a critical role in determining the permeability of the fractured reservoir. We found similar results in the calculation of equivalent permeability through topology analysis of the real cases. Larger branches of fracture systems tend to have more connectivity with other fracture planes, but I-node fractures tend to form short planes.

The shape factor was first introduced by Barenblatt et al. [84]. In this approach, the matrix and the fractures are divided into two continuum systems: The fractures are considered the main global flow paths (fractures have high permeability and low storage capacity), while the rock matrix is the main fluid storage sink (the matrix blocks present high storage and low permeability). The two systems are connected, and they interact either directly or indirectly with the global connection neighboring fractures. The concept of double porosity was extended and applied to reservoir engineering by Warren and Root [46], mainly for pressure test analysis. They developed a method based on the general premise that fractures are evenly spaced and allow for variations in fracture width to satisfy real anisotropy conditions in a naturally fractured reservoir. The fracture-matrix transfer term, governed by shape factor, cannot be represented by pseudo-steady-state. Rather, it involves transient periods of variable magnitude depending on the considered flow case. Furthermore, the shape factor is not just a value. The deficient concept of the pseudo-state transfer assumption explains why a certain controversy developed around the definition of the shape factor [85]. The fracture-matrix transfer problem involves a time convolution between fracture flow conditions. The actual flow occurs within the matrix medium.

Certain works have developed transient methodologies to obtain the shape factor, regarding the time effect, and a real representation of the fracture matrix interaction regarding the production or injection processes in the reservoir [23,47,51,52,53,86,87,88,89,90]. Rostami [91] demonstrated that the shape factor is a function of time, and its value could be different for regular and irregular shapes.

The analysis performed in this study considered a novel approach to obtain the shape factor as a function of time and compared the equivalent permeability computed through numerical simulation with analytical methods corrected with topological analysis in bidimensional artificial and real scenarios. The most relevant result demonstrated the influence of the I-node topology feature on the flow in fractured media. Fracture networks dominated by I-node types hinder the transfer effect between the fractures and the rock matrix, and the impact in the equivalent permeability can be relevant, despite the fracture density. Thus, the arrangement of fractures that influences the topological nodes and the hydraulic connectivity can provide a simplified and fastened estimate of fluid control parameters (permeability and shape factor) for numerical simulations of fractured reservoirs.

6. Conclusions

Numerical simulations were successfully performed to obtain the effective permeability regarding artificial and real cases with dual-porosity in fractured media built with FEM models. The results were compared with data obtained through analytical methods, corrected with topological analysis of fracture systems (I-nodes, Y-nodes, and X-nodes).

Results showed the impact of connection between the fractures and the behavior of properties directly related to flow, such as permeability. The analysis demonstrated the relationship between the fracture topology and shape factor. Fractured media dominated by I-node terminations result in smaller shape factors than systems composed by other fracture node types. This finding is directly related to the lack of connectivity.

Simulations performed on the scenarios showed that in fractured systems with fewer X- and Y-nodes than I-nodes, the equation used to define the fracture connectivity f will return negative results. This finding prevented the further extraction of properties such as equivalent permeability or the shape factor regarding the hydraulic conductivity f of the medium.

The results presented demonstrated the direct relationship of the shape factor with the geometry of the fractures and the fracture intensity aspect and evidenced the impact of these variables on the shape factor.

The comparative analysis demonstrated that the analytical methods are satisfactory to estimate equivalent permeability of idealized fractured systems. These represent artificial cases in terms of the fracture density and the connectivity factor. However, the test with real scenarios, a DFN from a pre-salt reservoir, demonstrated that the analytical methods are not efficient to estimate the equivalent permeability of the fractured medium, including with the corrected analytical method by using topological analysis.

The study reinforces the contribution of the fracture topology theory in obtaining input parameters for double porosity and double permeability models such as permeability and shape factor.

As a recommendation for the work evolution, one can consider the variable aperture and cement filling due to diagenetic effects in fractures, which is crucial to reproduce realistic models. The adoption of a coupled hydro-mechanical approach should be considered in further developments to establish the effect of fractured media (matrix rock and fractures) and stress-strain in the upscaled process need to be analyzed.

Author Contributions

Conceptualization, I.F.G., L.J.d.N.G., T.S.d.M. and L.L.A.; methodology, I.F.G., L.J.d.N.G., L.B., B.M., L.C.P., T.S.d.M., J.A.B. and L.L.A.; software, L.J.d.N.G., L.B., B.M. and L.L.A.; formal analysis, I.F.G., J.A.B., L.L.A. and T.S.d.M.; investigation, L.L.A., I.F.G., J.A.B., T.S.d.M., L.J.d.N.G. and B.M.; writing—original draft preparation, L.L.A., I.F.G., L.J.d.N.G., J.A.B. and T.S.d.M.; writing—review and editing, L.L.A., J.A.B. and I.F.G.; supervision, I.F.G., L.J.d.N.G. and J.A.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

The authors thank the financial support provided by Petrobras, Energi Simulation, and CNPq-Brazil.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Classification of node types. X type: fractures are crossed, I type: fracture start and end, Y type: fracture abutment with another fracture. Modified from Sanderson and Nixon [38].

Figure 2. (a) Finite element crossed by a discontinuity band. (b) The normal and tangential components of fluid flow in fracture.

Figure 3. Synthetic scenarios of natural fracture systems (a) I-node; (b) X-node; (c) Y-node; and (d) combination of XYI-nodes.

Figure 4. Synthetic models with I-, Y-, and X-node combinations and boundary conditions were used in the numerical simulation. (a) The horizontal fluid flow and (b) vertical flow.

Figure 5. (a) Geocellular model with a zone filled with fracture planes (DFN) developed for a pre-salt reservoir [32]. (b) Location of a grid-cell extracted with a sub-sampling of DFN model. (c) DFN sample used for the extraction of 2D sections used in this study.

Figure 6. Bidimensional sections extracted from the 3D model cell (slices). (a–c) sections 1, 2, and 3, with 54 fractures, 22 fractures, and 75 fractures, respectively.

Figure 7. Equivalent permeability output obtained through the methods tested regarding the simulation of horizontal flow in the synthetic cases (Figure 3).

Figure 8. Equivalent permeability output obtained through the methods tested regarding the simulation of vertical flow in the synthetic cases (Figure 3).

Figure 9. Computation of shape factor behavior with time for the synthetic scenarios.

Figure 10. Computation of shape factor behavior with time for the synthetic scenarios.

Figure 11. Cross-correlation between computed equivalent permeability and computed shape factor of synthetic cases.

Figure 12. Equivalent permeability output obtained through the methods tested regarding the simulation of horizontal flow in the three real cases from the pre-salt reservoir model (Figure 6).

Figure 13. Equivalent permeability output obtained through the methods tested regarding the simulation of vertical flow in the three real cases from the pre-salt reservoir model (Figure 6).

Figure 14. Cross-correlation plot between computed equivalent permeability and the hydraulic connectivity factor found through the topological analysis regarding the horizontal flow for the three bidimensional sections modeled after the pre-salt reservoir.

Figure 15. Cross-correlation plot between computed equivalent permeability and the hydraulic connectivity factor found through the topological analysis regarding the vertical flow for the three 2D sections modeled after the pre-salt reservoir.

Figure 16. Graph showing prediction of shape factor behavior with time for the three real sections after the pre-salt reservoir model.

Figure 17. (a) Comparison between fracture density, P22 (relationship of fracture area over the matrix area), (b) Computed shape factor for the three 2D sections extracted from the pre-salt reservoir model.

Figure 18. Cross-correlation between computed shape factor and hydraulic connectivity f for the three 2D sections of the pre-salt reservoir model.

Figure 19. Cross-correlation between equivalent permeability and the shape factor regarding the X- and Y-node types fractures for the three real scenarios of the pre-salt reservoir model.

Figure 20. (a) Ternary diagram with the relationship between the topological patterns I, Y, X, used to classify the DFNs according to their connectivity and (b) Comparison of hydraulic connectivity between sections.

Table 1

Quantity of fracture node type connections of synthetic scenarios 1, 2, 3, and 4.

Section	Node I	Node X	Node Y	Node XYI
Fractures	115	9	26	24
Pattern I	230	-	-	11
Pattern X	-	13	-	14
Pattern Y	-	-	34	15

Table 2

Quantity of fracture node type connections of sections 1, 2, and 3 extracted from the Brazilian pre-salt reservoir model cell (Figure 5).

	Section 1	Section 2	Section 3
Fractures	54	22	75
Pattern I	43	29	77
Pattern X	79	12	123
Pattern Y	12	3	23

Table 3

Summary of synthetic cases characteristics and the computed equivalent permeability and shape factor related to their hydraulic connectivity f.

Section	Node I	Node X	Node Y	Node XYI
Fractures	115	9	26	24
Pattern I	230	-	-	11
Pattern X	-	13	-	14
Pattern Y	-	-	34	15
Hydraulic Connectivity f	0.0	0.81	0.81	0.48
Matrix k (mD)	0.1	0.1	0.1	0.1
Keq x–Numerical (mD)	0.0107	0.1151	0.1153	0.1154
Keq x–Oda’s (mD)	0.0506	0.1390	0.1490	0.1406
Keq x–f correction (mD)	0.0	0.1126	0.1149	0.0670
Shape Factor	6.54	25.6	14.5	4.64

Table 4

Summary of real cases characteristics and the computed equivalent permeability and shape factor related to their hydraulic connectivity f.

	Section 1	Section 2	Section 3
Fractures	54	22	75
Pattern I	43	29	77
Pattern X	79	12	123
Pattern Y	12	3	23
Hydraulic Connectivity f	0.48	0.0	0.44
Matrix k (mD)	343.28	343.28	343.28
Keq x–Numerical (mD)	369.44	375.21	410.45
Keq x–Oda’s (mD)	407.37	409.56	464.52
Keq x–f correction (mD)	195.51	0.0	205.27
Keq y–Numerical (mD)	544.67	385.64	586.92
Keq y–Oda’s (mD)	705.02	488.73	913.90
Keq y–f correction (mD)	338.36	0.0	403.85
Shape Factor	0.1964	0.2766	0.0689