Evaluating Flow in Fractal-Fracture Networks: Effect of Variable Aperture

While fractal models are often employed for describing the geometry of fracture networks, a constant aperture is mostly assigned to all the fractures when such models are flow simulated. In nature however, almost all fracture networks exhibit variable aperture values and it is this fracture aperture that controls the conductivity of individual fractures as described by the well-known cubic-law. It would therefore be of practical interest to investigate flow patterns in a fractal-fracture network where the apertures scale in accordance to their position in the hierarchy of the fractal. A set of synthetic fractal-fracture networks and two well-connected natural fracture maps that belong to the same fractal system are used for this purpose. A set of dominant sub-networks are generated from a given fractal-fracture map by systematically removing the smaller fracture segments with narrow apertures. The connectivity values of the fractal-fracture networks and their respective dominant sub-networks are then computed. Although a large number of fractures with smaller aperture are eliminated, no significant decrease is seen in the connectivity of the dominant sub-networks. A streamline simulator based on Darcy's law is used for flow simulating the fracture networks, which are conceptualized as two-dimensional fracture continuum models. A single high porosity value is assigned to all the fractures. The permeability assigned to fractures within the continuum model is based on their aperture values and there is nearly no matrix porosity and permeability. The recovery profiles and time-of-flight plots for each network and its dominant sub-networks at different time steps are compared. The results from both the synthetic networks and the natural data show that there is no significant decrease in fluid recovery in the dominant sub-networks compared to their respective parent fractal-fracture networks. It may therefore be concluded that in the case of such hierarchical fractal-fracture systems with scaled aperture, the smaller fractures do not significantly contribute to connectivity or fluid flow. In terms of decision making, this result will aid geoscientists and engineers in identifying only those fractures that ultimately matter in evaluating the flow recovery, thus building models that are computationally less expensive while being geologically realistic.

1 Introduction

Fluid flow and transport characteristics in fractured media are controlled by connectivity and conductivity of fracture networks. The effect of connectivity on flow behavior of fracture networks is well studied by many researchers (Robinson, 1983; Berkowitz and Balberg, 1993; Berkowitz, 1995; de Dreuzy et al., 2002; Sahu and Roy, 2020). However, when it comes to simulating flow in such networks, an average value of aperture is assigned in most cases (Min et al., 2004; Hardebol et al., 2015; Chen, 2020) even though it has been recognized that such assumption of constant fracture aperture is an oversimplification (Odling and Roden, 1997). Aperture is the void space between two surfaces of a fracture. They vary in size from micrometres to centimetres as they result not only from mechanical misfit of fracture walls but also from chemical dissolution, and overburden pressure (Bonnet et al., 2001). Hence, the characterization of aperture distribution is important as it controls the rock mass hydro-mechanical properties of fracture networks (Kaulatilake et al., 2006). Neretnieks et al. (1982) experimentally illustrated that flow channelling occurs in individual fractures and mostly through highly preferential pathways in fracture networks. Cacas et al. (1990a, b) showed that broad distribution of fracture conductivities (apertures) is the main cause of high degree of flow channelling. Berkowitz (2002) pointed out that even a well-connected fracture network can display sparse preferential flow paths if the distribution of fracture aperture is significantly wide. Klepikova et al. (2014) showed that inversion of borehole temperature profiles under predefined pumping conditions can be applied for the detection of main permeable fractures. Although this approach of using temperature tomography does not identify all flowing fractures, it allows the detection of the most transmissive ones. A study on fracture insertion and deletion based on DFN transdimensional inversion using reversible jump Markov Chain Monte Carlo (rjMCMC) algorithm was carried out by Somogyvari et al. (2017). One interesting conclusion of this research is that no breakthrough is observed in fractures with smaller apertures within the simulated time. However, there is limited information in the literature that help researchers in realistically assigning a distribution of apertures in fracture network models.

In many cases, such fracture aperture distribution is observed to follow a power law in naturally fractured reservoirs (Bour and Davy, 1997) and hence, the occurrence of smaller fractures with narrow aperture is more frequent than fractures of longer length with broader aperture (Somogyvari et al., 2017). The fluid flow behavior in each individual fracture element is controlled by the hydraulic aperture and length of the fractures (Bour and Davy, 1997; de Dreuzy et al., 2002; Darcel et al., 2003; Baghbanan and Jing, 2007). Therefore, the characterization of aperture distribution in fracture networks and its effect on flow behaviour is crucial in identifying the dominant fractures that serve as preferential pathways in fracture networks.

This research investigates how aperture distribution affects the connectivity and fluid flow in fractal-fracture networks and whether fractures with smaller apertures contribute significantly in terms of the overall flow recovery. One of the objectives of this study is to aid modellers in their decision-making process on how much of information (data) should be incorporated in building fracture models that are geologically realistic, i.e., the smallest apertures that can be eliminated from models without significantly affecting the flow while being computationally less expensive. In order to meet this objective, dominant sub-networks have been identified that are generated by the elimination of smaller fractures from a parent fractal-fracture network but nonetheless, have very similar fluid recovery to that of the latter. Fractures are assumed to have clean openings, that is, those that have the entire fracture width available for flow. In other words, we consider the hydraulic aperture (Snow, 1969; Neuzil and Tracy, 1981; Barton et al., 1985).

The connectivity and flow behavior of a set of synthetic fractal-fracture networks, both deterministic and random (Roy et al., 2007), and their dominant sub-networks, where apertures scale in accordance to their position in the hierarchy of a fractal model, are studied. A dominant sub-network is obtained by systemically eliminating fractures with smaller aperture from the parent network (Gong and Rossen, 2017). In order to simulate flow, a fracture continuum (FC) model is considered which is similar to those used by Langevin (2003), Neuman (2005), Svensson (2001a, b) and Tsang et al. (1996). Cells in the FC model are assigned porosity and permeability values depending on whether they represent fracture or matrix. A high porosity value is assigned to all the cells that represent fractures. Permeability values of “fracture cells” are computed based on the cubic law (aperture) and there is nearly no matrix porosity or permeability. Trace3D, a Darcy based streamline simulator (Datta-Gupta and King, 2007) is used for modelling fluid flow in the FC models. The concept thus developed, is then applied to a natural fracture map (Odling, 1997) and its dominant sub-network to demonstrate that in case of real data, fractures with smaller apertures may be excluded from the flow model without significantly affecting recovery. The dominant sub-network in this case has been generated from a lower resolution map of original fracture map which captured only the thick and longer fractures.

2 Method Development

2.1 Generation of fractal-fracture networks and their dominant sub-networks

Fractal-fracture networks can serve as useful models for quantitative evaluation of heterogeneity and complexity of reservoir rocks. These models can be used to study properties like scaling of fractures, aperture distribution and fracture network geometry that affect the fluid flow behavior in a fracture system (Berkowitz and Hadad, 1997). In this study, fractal-fracture networks are modeled as deterministic and random Sierpinski lattices composed of self-similar segments in two dimensions. These are generated using three parameters: the scale factor, b, the iteration level, i, and the initial number of un-fractured blocks, n (here b=2, i=5 and n=1) as described in Roy et al. (2007). In order to investigate the effect of aperture on fluid flow in such networks, rather than using a constant value, the fracture apertures are assigned a distribution such that they scale in accordance with their hierarchy in the network as described in Eq. (1). It may be noted that there exists a positive correlation between the length of the fractures and their aperture values in this model.

$$\begin{array}{}\text{(1)}& w_i=b^{i\left(\max \right)-i}\end{array}$$

Here, b is the scale factor, w_i is the aperture of a fracture segment that was generated at the iteration, i, and, i_max is the maximum number of iterations used for generating the entire fractal-fracture network. In this case, b=2 and $i_{max}=\mathrm{5}$.

A set of dominant sub-networks (Fig. 1) are generated from a given parent fractal-fracture network by systematically eliminating fractures with smaller apertures from the latter (Gong and Rossen, 2017). This is done in order to test whether fractures with such smaller apertures significantly contribute to the connectivity and flow in such networks. It helps in understanding whether the sub-networks (left after removing the smaller apertures) carry most of the injected fluids. Further, it also addresses the importance of these sub-networks in characterization of flow in fractured reservoirs.

Figure 1 (a–c) Deterministic fractal-fracture network (a) and its dominant sub-network-I (b) and dominant sub-network-II (c). (d–f) Random fractal-fracture network (d) and its dominant sub-network-I (e) and dominant sub-network-II (f).

As the dominant sub-networks are obtained by removing fractures with smaller aperture, it is important to keep count of the numbers of fractures eliminated at each step. A MATLAB toolbox, FracPaQ (Healy et al., 2017) is used for this and the results are shown in Fig. 2. Compared to that of the parent fractal-fracture network, the number of fractures decreases by 75 % in dominant sub-network-I and by 85 % in dominant sub-network-II.

Figure 2 Number of fracture segments vs. aperture for a fractal-fracture network and its dominant sub-networks-I and -II.

In Fig. 2, it is seen that the parent fractal-fracture network (red bars) consists of apertures of all values, i.e., 100 to 3200 µm. Dominant sub-network-I (green bars) is generated by removing the 100 µm apertures from the parent network and has apertures ranging from 200 to 3200 µm. Dominant sub-network-II (purple bars) is generated by removing both the 100 and 200 µm apertures. It therefore, has apertures ranging from 400 to 3200 µm. Figure 2 is applicable to both deterministic and random fractal-fracture networks as they have the same length and aperture distributions but differ only in the spatial distribution of the fracture segments.

2.2 Connectivity of fractal-fracture networks and their dominant sub-networks

The connectivity of a fracture network is an important factor in determining flow and transport characteristics (Robinson, 1983; Berkowitz and Balberg, 1993; Berkowitz, 1995; de Dreuzy et al., 2002; Sahu and Roy, 2020). Topological graph theory uses the relative abundance of three types of nodes in a fracture network to define the connectivity (Sanderson and Nixon, 2015; Sanderson et al., 2019). These nodes are: X (cross-cutting), Y (abutting) and I (isolated) as shown in Fig. 3. The relative abundance of these different types of nodes in fractal-fracture networks and their respective dominant sub-networks are estimated by FracPaQ, a MATLAB tool box for quantification of fracture patterns from 2D images (Healy et al., 2017).

Figure 3 Schematic fracture network showing three different types of nodes.

Manzocchi (2002) defines the topological connectivity of fracture networks in terms of a single parameter, n as in Eq. (2):

$$\begin{array}{}\text{(2)}& n=\left[\mathrm{4}\right(\mathrm{1}-P_{\mathrm{I}})/(\mathrm{1}-P_{\mathrm{X}}\left)\right]\end{array}$$

Here, P_I, P_Y and P_X are the proportions of I (isolated) nodes, Y (abutting) nodes and X (crosscutting) nodes respectively. The sum of P_I, P_Y and P_X is always equal to 1.

The connectivity values of the deterministic and random fractal-fracture networks and, their respective dominant sub-networks, are shown in Table 1. It is no surprise that the connectivity of the dominant sub-networks decreases because these are generated by eliminating a large number of smaller aperture fractures from the parent network. However, in the case of deterministic network, this decrease is only about 12 % in dominant sub-network-I and 25 % in dominant sub-network-II. For the random fractal-fracture network, it is decreased by 10 % in dominant sub-network-I and 23 % in dominant sub-network-II. It may be noted that the connectivity does not decrease by more than 25 %, even after elimination of a large number of smaller fractures. This is because, elimination of the smaller fractures from a parent fractal-fracture network results in a decrease in all three types of nodes by a similar proportion to what is present in the parent network, thus keeping the proportion of nodes almost same in the resultant, dominant sub-network.

Table 1 Connectivity values of different fracture networks and their dominant sub-networks.

2.3 Flow simulation in fractal-fracture networks and their dominant sub-networks

Simulating flow through fracture networks with variable aperture is computationally demanding, especially when it comes to using DFN models (Schwartz and Smith, 1988; Parney and Smith, 1995; Benke and Painter, 2003; Painter and Cvetkovic, 2005; Leung et al., 2012). Fracture continuum (FC) models on the other hand, are computationally efficient while preserving the network details and fracture properties. FC models have been widely used for numerical flow simulations in fractured rocks (Langevin, 2003; Neuman, 2005; Svensson, 2001b; Tsang et al., 1996; McKenna and Reeves, 2006; Botros et al., 2008; Reeves et al., 2008).

In FC models, discrete fractures are modelled as permeability structures on a grid as shown in Fig. 4. In this study, the synthetic fractal-fracture patterns and their dominant sub-networks are modelled as grids of $\mathrm{500}\times \mathrm{500}\times \mathrm{1}$ cells, such that each pixel of the network's image (raster data) is represented by a cell. The details of the fracture network are preserved by the use of a very small cell size, while a high permeability contrast between fractures and matrix can restrict the flow only within the former.

Figure 4 Illustration of finite difference representation of a hypothetical rock fracture network (modified from Reeves et al., 2008).

Following the footsteps of Sarkar et al. (2004) the permeability, K_f, of cells representing fractures, is assigned based on the well-known “cubic law” (Snow, 1969) as in Eq. (3):

$$\begin{array}{}\text{(3)}& K_f=w^{\mathrm{2}}/\mathrm{12}\end{array}$$

Here, w is the hydraulic aperture. The values of fracture permeability calculated for given fracture apertures are shown in Table 2. Cells that represent matrix (white cells in Fig. 4) are assigned negligible permeability and porosity such that the flow mainly occurs in the fractures.

Table 2 Fracture aperture and permeability values.

The fractal-fracture networks and their respective dominant sub-networks are flow simulated using Trace3D (Dutta-Gupta and King, 2007). It is a streamline simulator which represents fluid flow in a reservoir by implementing Darcy's law and mass conservation equations. In this simulator, approximate changes in pressure and saturation distribution with time are obtained by solving pressure and saturation/conservation equations (Al-Najem et al., 2012). Such an approach is well-suited for the analysis of fluid flow characteristics with recovery curves and time of flights (TOFs) as outputs. The IMPES (implicit pressure explicit saturation) method allows the pressure equations to be solved implicitly, i.e., pressures at all grid nodes are computed simultaneously by solving a system of linear equations at every time step. On the other hand, IMPES computes saturations explicitly, i.e., one grid block at a time without any matrix solution. The solution of the pressure and saturation generates the fluid recovery values and also maps the advancement of fluid front at different time steps represented by TOFs. The details of all the governing equations used in Trace3D can be found in Datta-Gupta and King (2007).

An injection and a production well are placed at diagonally opposite corners of the model grid (i.e., 1×1 and 500×500) to obtain a total areal sweep of the fractured rock volume. Flow simulation is performed at a reservoir pressure of 2480 psia with the assumption of constant boundary reservoir conditions. The fluid injection rate is maintained at 500 bbl/day for a simulation time period of 700 d. The entire fluid volume is considered to be stored in the fractures and therefore, a porosity value of 95 % is assigned to cells that represent fractures. The permeability of “fracture cells” is assigned depending on their apertures as discussed earlier and shown in Table 2. The porosity and permeability of “matrix cells” are assigned as 0.05 % and $\mathrm{0.98}\times {\mathrm{10}}^{-\mathrm{16}}$ m² respectively, such that flow through the matrix is negligible.

3 Results and Discussion: Synthetic Fractal-Fracture Networks

Fluid recovery can serve as a parametric estimation of flow characteristics in any fracture network. Figure 5 shows the recovery curves for the deterministic and random fractal-fracture networks and their respective dominant sub-networks. For the deterministic network (solid lines in Fig. 5) it can be seen that the dominant sub-networks are almost indistinguishable from their parent network in terms of recovery. A similar behavior is seen for the random network (dashed lines in Fig. 5), although some differences exist in the dominant sub-networks in this case. This indicates that even though the parent fractal-fracture networks and their dominant sub-networks have very different number of fractures (Fig. 2), they have almost similar fluid flow characteristics. It is also seen from the figure that the response from dominant sub-network-I of the deterministic pattern actually coincides with that of the parent network of the random pattern. As because the random and deterministic patterns have the same length distribution it may be conjectured that the parent fractal-fracture network is very similar to the sub-networks in terms of recovery.

The overall recovery values thus obtained for all fracture networks studied are shown in Table 3. Fluid recovery for dominant sub-network-I is about 96 % and 93 % to that of the parent network for the deterministic and random patterns respectively. In case of dominant sub-network-II, these values are about 85 % and 83 %. Hence, it may be hypothesized that while modelling such fracture networks, eliminating a large number fractures with smaller apertures will not significantly affect the fluid flow.

Figure 5 Recovery curves of deterministic and random fractal-fracture networks along with their respective dominant sub-networks.

Table 3 Overall fluid recovery values of different fracture networks and their dominant sub-networks.

Time of flight (TOF) is an important parameter for visualizing the fluid velocity at a given location within the flow domain. The TOF can be defined as a spatial coordinate representing the distance along streamline to where the reservoir will be contacted at that instance of time. It is worth mentioning that along preferential flow paths, the TOFs show higher resolution as the streamline's tends to cluster along the high permeability streaks. The TOF's depicting fluid-front through the deterministic and random fractal-fracture network along with their dominant sub-networks at different time steps (500, 1500 and 3000) are shown in Fig. 6. It is observed that flow patterns are almost similar for the parent fractal-fracture network and their dominant sub-networks. The TOFs clearly show flow channelling through fractures with larger apertures of the dominant sub-networks in both deterministic and random fractal-fracture networks. It can be visually seen from Fig. 6 that these fractures serve as preferential pathways for conduit of major amount of fluid in the reservoir and is in agreement with the experimental analysis of Neretnieks et al. (1982). The time of flight for fractures with larger apertures are much smaller than those with smaller apertures hence, fluid movements are faster in the former. This results in almost similar recovery in the original fractal-fracture network and its dominant sub-networks which contain only the larger fractures.

Figure 6 (a) TOFs at different time periods of deterministic fractal-fracture network (top), its dominant sub-network-I (centre) and -II (bottom). (b) TOFs at different time periods of random fractal-fracture network (top), its dominant sub-network-I (centre) and -II (bottom).

4 Application to Natural Fracture Networks

A subset of natural fracture networks is chosen from a set of seven nested fractal-fracture maps gathered from the Devonian Sandstones of Hornelen basin, Norway (Odling, 1997). These maps have been studied by many researchers in relation to the characterization of various statistical and geometrical properties of fracture networks in the last three decades (Bonnet et al., 2001; Bour et al., 2002; Odling and Roden, 1997; Roy et al., 2007, 2010). Hence, they may serve as good outcrop analogs for fractured reservoirs. They cover areas from 18 to 720 m across at a range of resolutions forming a nested series. Each map represents a “window” or sub-sample of the fracture trace population with a lower bound controlled by the resolution and an upper bound controlled by the size of the region observed (Bour et al., 2002; Odling, 1997; Roy et al., 2007). The scale, resolution and the height of observation of these maps along with the mapping technique are described in Odling (1997). We chose map 1 and map 2 from the set such that the latter contains the former but, at a lower resolution as shown in Fig. 7. The aperture of the natural fracture map is scaled considering the trace-length distribution of individual maps which follows a power law (Odling, 1997).

Figure 7 Illustration of mapping technique of natural fracture map and its dominant sub-network map (modified from Odling, 1997).

A dominant sub-network of map-1 is generated by considering a smaller section from map-2 that represents map-1 (at a lower resolution) and rotating it to match its orientation with the latter. Thus, the dominant sub-network of map-1 retains only the longer fractures (with large apertures) as the shorter fractures (with smaller apertures) are “lost” due to decrease in resolution. The count of fractures in both the maps is estimated using FracPaQ and the dominant sub-network contains only about 15 % of the fracture segments from the original parent fracture network. The natural fracture network (map-1) and its sub-dominant network are shown in Fig. 8.

Figure 8 Natural fracture map (left) and its dominant sub-network (right) (modified from Odling, 1997).

The proportion of X, Y and I nodes present in the natural fracture map and its dominant sub-network are estimated using FracPaQ and their connectivity values are calculated from Eq. (2). The connectivity of the dominant sub-network is observed to decrease by 15 % of that of its parent fracture map as shown in Table 2. The natural fracture map and its dominant sub-network are modeled as grids of $\mathrm{1042}\times \mathrm{1042}\times \mathrm{1}$ cells invoking the fracture continuum (FC) concept as before. Similar to our approach in Sect. 2.3, the permeability of cells that represent fractures are calculated from the cubic law, and a single high porosity of 95 % is assigned to all cells representing fractures. As with the models, matrix porosity and permeability values are assigned as 0.05 % and $\mathrm{0.98}\times {\mathrm{10}}^{-\mathrm{16}}$ m² respectively, such that flow through the matrix is negligible. The other inputs used in the simulator are also kept the same as mentioned in Sect. 2.3. The fracture maps are flow simulated using Trace3D, with an injection and a production well placed at diagonally opposite corners of the model grid (i.e., 1×1 and 1042×1042). The recovery plots of the natural fracture map and its dominant sub-network are shown in Fig. 9. Their TOFs that represent fluid movement at different time steps (500, 1500 and 3000) are shown in Fig. 10.

Figure 9 Recovery profile of scaled natural fracture map and its dominant sub-network.

Figure 10 TOFs at different time steps (500, 1500 and 3000) of natural fracture map (a–c) and its dominant sub-network (d–f).

The overall recovery values of each of the two maps are shown in Table 3. The fluid recovery value for the dominant sub-network is almost 88 % of that of the original map although the former retains only 15 % of the fractures from that of the latter. This clearly indicates that fractures with smaller aperture have a rather negligible contribution to the overall fluid flow. The TOFs for the natural fracture map and its dominant sub-network show almost no difference at least up to 1500 time steps. Some minor differences are observed between the two networks at 3000 time steps. This shows that the smaller fractures present in the original map (and those that are captured only in high resolution data) actually do not significantly contribute to the flow. From a modelling perspective, this observation is very helpful because geoscientists need not include data on smaller fractures in their models as they do not significantly affect the connectivity and overall fluid recovery.

5 Concluding Remarks

The influence of aperture distribution on connectivity and fluid flow in fractal-fracture networks has been explored by studying a set of model networks, a natural map and their respective dominant sub-networks. Although the dominant sub-networks were generated by removing a large number of smaller aperture fractures from the parent networks, the connectivity was observed to decrease only by about 10 %–25 %. The networks were then flow simulated in Trace3D invoking the FC model. It was seen that in networks with a broad distribution of aperture, most of the fractures with smaller apertures can be eliminated without any significant reduction of overall fluid recovery. This is because dominant subnetworks containing only about 15 % of the total fracture segments from their respective parent networks yielded recovery values nearly 85 % of that of the latter.

The same strategy was applied to a natural fracture map where smaller fractures were removed by comparing it to its lower-resolution version such that only the long and wide fractures that appeared in the latter were retained. Here, the connectivity and fluid recovery of the dominant sub-network map was found to be around 84 % and 88 % of the parent network respectively. Thus, the results from this experiment proved to be very similar to what was obtained from the study of model fractal-fracture networks. It may therefore be concluded that even in a well-connected fracture network, fluid flow may be unequally distributed based on the aperture distribution and most of the network can be excluded, 85% in this particular case, without significantly affecting its fluid recovery.

In terms of decision making, this result will aid geoscientists in eliminating smaller fractures from the data that do not contribute significantly to the flow thus, building fracture models that are computationally less expensive. The implications of this finding for characterization of fractured rock mass can prove useful in primary as well as enhanced methods of fluid recovery, where the dominant sub-network fractures are the one that carry most of the produced and injected fluids as the larger aperture fractures provide preferential pathways for fluid flow.