Shape factor for dual-permeability fractured reservoir simulation: Effect of non-uniform flow in 2D fracture network

The flow properties of naturally fractured reservoirs are dominated by flow through the fractures. In a previous study we showed that even a well-connected fracture network behaves like a much sparser network when the aperture distribution is broad enough: i.e., most fractures can be eliminated while leaving a sub-network with virtually the same permeability as the original fracture network. In this study, we focus on the influence of eliminating unimportant fractures which carry little flow on the inferred characteristic matrix-block size. We model a two-dimensional fractured reservoir in which the fractures are well-connected. The fractures obey a power-law length distribution, as observed in natural fracture networks. For the aperture distribution, because information from the subsurface is limited, we test a number of cases: log-normal distributions (from narrow to broad), and power-law distributions (from narrow to broad). The matrix blocks in fractured reservoirs are of varying sizes and shapes; we adopt the characteristic radius and the characteristic length to represent the characteristic matrix-block size. We show how the characteristic matrix-block sizes increase from the original fracture network to the dominant sub-network. This suggests that the matrix-block size, or shape factor, used in dual-porosity/dual-permeability waterflood or enhanced oil recovery (EOR) simulations or in homogenization should be based not on the entire fracture population but on the sub-network that carries almost all of the injected fluid (water or EOR agent).


Introduction
Naturally fractured reservoirs contain a significant amount of hydrocarbon reserves worldwide [1], However, the oil recovery from these reservoirs has been rather low. The low level of oil recovery indicates that more accurate reservoir characterisation and flow simulation is needed.
Reservoir simulation is one of the most practical methods of studying flow problems in porous media. For fractured reservoir simulation, the dual-porosity/dual-permeability concept and the discrete fracture model are two typical methods [2]. In the dualporosity/dual-permeability approach, the fracture and matrix systems are treated as separate domains; interconnected fractures serve as fluid flow paths between injection and production wells, while the matrix acts only as fluid storage, and these two domains are connected with an exchange term [3][4][5]. In a dual-permeability model, fluid flow can also take place between matrix grid blocks, unlike from the dual-porosity model [6,7]. In order to simulate the realistic fracture geometry and account explicitly for the effect of individual fractures on fluid flow, discrete-fracture models have been developed [8][9][10][11][12][13][14]. Compared to the dual-porosity/dualpermeability models, discrete-fracture models represent a fracture network more explicitly and make the simulation more realistic. But discrete fracture models are typically difficult to solve numerically. Thus, although dual-porosity/dual-permeability models are much simplified characterizations of naturally fractured reservoirs, they are still the most widely used methods for field-scale fractured-reservoir simulation, as they address the dual-porosity nature of fractured reservoirs and are computationally cheaper. To generate a dual-porosity/dual-permeability model, it is necessary to define average properties for each grid cell, such as porosity, permeability and matrix-fracture interaction parameters (typical fracture spacing or shape factor) [15]. Therefore, the discrete fracture network considered to generate the dual-porosity model parameters is crucial. Using homogenization, one can treat matrix-fracture exchange more accurately than in dual porosity/dual permeability simulations [16], but, again, one needs a characteristic matrix-block size.
As we presented in a previous study [17], even in a wellconnected fracture network, there is a dominant sub-network which carries almost all the flow, but it is much more sparse than the original network (Fig. 1). The flow-path length of the dominant sub-network can be as little as roughly 30% of that of the corresponding original fracture network in the most extreme case. This suggests that in secondary production, the water injected flows mainly along a small portion of the fracture network. In contrast, in primary production even relatively small fractures can be an efficient path for oil to flow to a production well.
This paper is organized as follows: we first introduce the numerical model and research process. Next, we analyse the sizes of the matrix blocks formed by the entire fracture network and the corresponding dominant sub-network. Finally, we point out the implications of this distinction for the dual-porosity/ dual-permeability reservoir simulation.

Models
Since this is a follow-up study to our previous research [17], the models used here are the same as the ones adopted before (Fig. 1). Here we only introduce the models briefly; more details can be found in the previous work. Fracture networks are generated in a 10 m Â 10 m Â 0.01 m region using the commercial fracturedreservoir simulator FracMan TM [18]. Two fracture sets which are nearly orthogonal to each other are assumed, with almost equal numbers of fractures in the two sets. Each fracture, with a rectangular shape, is located following the Enhanced Baecher Model and is perpendicular to the plane which follows the flow direction and penetrates the top and bottom boundaries of the region. Because of the uncertainties in data and the influence of cut-offs in measurements, fracture-trace lengths have been described by exponential, log-normal or power-law distributions in previous studies [19][20][21]. Commonly, a power-law distribution is assumed by many researchers to be the correct model for fracture length [22][23][24][25], with exponent a ranging from 1.5 to 3.5. In this study, the fracture length follows the power-law distribution ðf ðxÞÞ: where a is the power-law exponent, x is the fracture length and x min the lower bound on x, which we take to be 0.2 m. In order to make sure that there are no extremely short or long fractures, and in particular that opposite sides of our region of interest cannot be connected by a single fracture, we choose a ¼ 2 and truncate the length distribution on the upper end at 6 m. Since even the smallest fracture length (0.2 m) is much larger than the thickness of the region of interest (0.01 m), the 3D model can be seen as a 2D fracture network.
For fracture apertures, two kinds of distribution which have been proposed in previous studies are adopted: power-law [26][27][28][29][30][31] and log-normal [32][33][34][35][36][37]. In each kind of distribution, to include the entire range of feasible cases (from narrow to broad aperture distribution), different parameter values (a for a power-law aperture distribution and r for a log-normal aperture distribution) are examined. The aperture is randomly assigned to each fracture.
The power-law distribution can also be defined as: If the power-law aperture distribution is described by Eq. (2), the studies cited above found that the value of the exponent a in nature is 1, 1.1, 1.8, 2.2, or 2.8. In this study, the power-law aperture distribution is defined by Eq. (1) as well as the fracture length distribution, where x stands for aperture instead of length. Eq. (1) differs from Eq. (2), in that it includes a minimum cut-off value, and a should be larger than 1. To include the entire range of feasible cases (from narrow to broad aperture distribution), here we examine a in the range from 1.001 to 6. In each case, the fracture aperture is limited to the interval between 0.01 mm and 10 mm. The aperture range can vary greatly in different formations; it also depends on the resolution and the size of the region studied. According to the field data collected from the Ship Rock volcanic plug in NW New Mexico [38] and Culpeper Quarry and Florence Lake [39], the aperture range [0.01 mm, 10 mm] adopted here is realistic, at least at some locations in natural. The log-normal distribution is specified by the following probability density function: where l and r are the mean and the standard deviation in the log-10 space. The truncated log-normal distribution has two additional parameters: a minimum and a maximum value of aperture. Field studies and hydraulic tests found values of r from 0.1 to 0.3, 0.23, and 0.47 [32,35,40]. To test the widest range of feasible values, we test values of r from 0.1 to 0.6. More details can be found in our previous study [17]. In order to focus on the influence of fracture aperture distributions on the dominant sub-network, except for the aperture distribution, all the other parameter distributions remain the same for all the cases tested in this study, including the fracture length, the orientation, etc. We assume that a fracture can be approximated as the slit between a pair of smooth, parallel plates; thus the aperture of each fracture is uniform. Steady state flow through the fractured region is considered (Fig. 1a). In this paper, we consider that fracture permeability is much greater than matrix permeability, which is common in fractured reservoirs [41,42]. The flow regimes of fractured formations can be defined by the fracture-matrix permeability ratio. Especially, if the ratio is greater than 10 5 -10 6 , fractures carry nearly all the flow [43,44]. Since we are interested in the non-uniform flow in well-connected fracture networks for characterizing flow in the fracture network, the matrix is further assumed to be impermeable, so that fluid flow can take place only in the fracture network, similar to the flow regime between fractures and permeable rock mass with the fracture-matrix permeability ratio is greater than 10 5 -10 6 . Assuming zero matrix permeability is merely a convenient assumption for calculating flow in fracture networks. For computing flow in discrete fracture networks, as in most numerical simulation methods, Darcy's Law for steady-state incompressible flow is employed, and mass is conserved at each intersection of fractures. In our models, we induce fluid flow from the left side to the right side by applying a constant difference in hydraulic head across the domain, while all the other boundaries are impermeable. Mafic TM , a companion program of FracMan TM , is employed to simulate flow in the fracture networks.

Methodology
As presented in the previous study, even in a well-connected fracture network, when the aperture distribution is broad enough, there is a dominant sub-network which can be a good approximation of the actual fracture network. This dominant sub-network is also strongly affected by the aperture distribution. The ''dominant sub-network" is defined as the sub-network obtained by eliminating a portion of fractures while retaining 90% of the originalnetwork equivalent permeability There is no special meaning for this 90% threshold, we choose it to illustrate the effect of nonuniform flow. In this study, our main interest lies in examining the change of the characteristic sizes of matrix blocks as more fractures are eliminated. We also check the influence of aperture distribution (exponent a in a power-law distribution and standard deviation r in a log-normal distribution) on the characteristic sizes of matrix blocks formed by the dominant sub-network. A Matlab program is employed to calculate the characteristic matrix block sizes. The equivalent matrix block size is employed to represent the average value of the resulting distribution of matrix block sizes.
The approach used to decide which portion of fractures to remove is as follows. Mafic TM subdivides the fractures into finite elements for the flow calculations. The flow velocity at the centre of each finite element and the value of flow velocity Â aperture (=Q nodal ) can be obtained. Based on this value, we compute the average value (Q) of all the elements in each fracture. Q is then used as the criterion to eliminate fractures: fractures are eliminated in order, starting from the one with the smallest value of Q to the one with the largest Q. That is to say, the fractures that conduct the least flow are eliminated first. After each step, we calculate the equivalent network permeability of the truncated network.
Because the generation of the fracture network is a random process, an infinite number of fracture networks can be generated with the same parameter values for the density, the orientation, the fracture length and the aperture distribution. In this study, for each set of parameter values, we generate one hundred realizations.

Characteristic matrix block sizes
The matrix blocks in fractured reservoirs can be of varying shapes and sizes. In order to quantitatively study the size and the recovery behaviour of the matrix blocks, several parameters have been proposed [45][46][47][48][49]. In this study, we adopt the characteristic radius and the characteristic length.
Zimmerman and Bodvarsson [50] argued that an irregularlyshaped matrix block can be modelled with reasonable accuracy as a spherical matrix block. The effective radius of the corresponding spherical block is defined as where V is the volume of the matrix block, and S is the surface area. They also proposed that in the early-time regime (shortly after a change in the boundary condition imposed at the block surfaces), a series of spherical blocks with variable radii can be modelled using uniform blocks with an equivalent radius given by: where n is the number of matrix blocks, i refers to one matrix block, r i is the effective radius of that matrix block, V i is the volume of that matrix block and V t is the total matrix volume. In this study, we adopt a similar idea to define the characteristic matrix-block radii of 2D fractured formations. The characteristic radius of each irregularly-shaped matrix block is defined as: where i refers to one matrix block, S mi is the area of the matrix block, and L mi is the perimeter of the two-dimensional matrix block. The equivalent matrix-block radius is then defined as: where n is the number of matrix blocks, i refers to one matrix block, S i is the area of that matrix block and S t is the total matrix area. The second parameter used in this study for representing the matrix-block size is the characteristic length, which was first proposed by Kazemi et al. [46] and followed by other researchers [47,51,52]. Kazemi et al. [46] proposed a shape factor (F s ) of a single matrix block for the imbibition process, which considered the effect of matrix-block shapes and boundary conditions: where i refers to one matrix block, V mi is the volume of that matrix block, j refers to one face of that matrix block which is open to imbibition, A j is the area of that face, d j is the distance from that face to the centre of the matrix block, and n is the total number of faces of the matrix block open to imbibition. This shape factor is claimed to be valid for anisotropic matrix blocks with irregular shapes [53]. The characteristic length of an irregular matrix block L c is then defined as In fractured reservoirs, the matrix block is the main location for storage of oil and it feeds the surrounding fractures; thus the bulk volume of the matrix is vital to the recovery rate. Therefore, we believe that a volume-weighted equivalent length is more reasonable than a number-based average value. An equivalent length for a series of matrix blocks is then represented as where i refers to one matrix block, V i is the volume of that matrix block, L ci is the characteristic length of that matrix block, and V t is the total bulk volume of the matrix blocks.
In this study, the equivalent radius (r e ) and the equivalent length (L e ) defined above are used to represent the average size of the matrix blocks in a fractured region. Only fractures that belong to the spanning cluster are included in calculating the matrix block sizes; the dead-ends are systematically removed from the network. Also, the matrix blocks containing impermeable boundaries (i.e. along the top and bottom of the region of interest) are not considered. Only the matrix blocks formed by the fractures which have fluid flow through them and the permeable left and right boundaries are taken into account.

Results
In this section, we present the results for the cases with powerlaw aperture distributions (from narrow to broad) and log-normal aperture distributions (from narrow to broad).  (Fig. 2d). The normalized equivalent permeability of the truncated fracture network (K b ) is shown in Fig. 3 in which K b =K o is the ratio between the equivalent permeability of the dominant sub-network (K b ) and the equivalent permeability of the original fracture network (K o ), while l b =l o represents the ratio between the length of the backbone of the truncated fracture network (l b ) and the total length of the original fracture network (l o ). The results show that the dominant sub-network (the sub-network retaining 90% of the original equivalent network permeability after eliminating a portion of fractures) is strongly affected by the aperture distribution: for the broadest aperture distribution cases (a = 1.001), the flowpath length (l b ) is roughly 35% of the total length of the original fracture network (l o ), while for the narrowest aperture distribution (a = 6), the ratio of l b =l o is approximately 0.6 (Fig. 3). Correspondingly, as presented in Fig. 2d, the equivalent matrix radius for the dominant sub-network is on average about twice that for the original fracture network when a = 1.001 (the red dashed line), while the ratio is approximately 1.5 for the cases with the narrowest aperture distribution with a = 6 (the red dotted line). Fig. 4 shows the characteristic radii (r c ) of all the individual matrix blocks formed by the entire fracture network and the dominant sub-networks (retaining 90% of original-network equivalent permeability) for one realization for each value of a (Fig. 4). The vertical axis (''CDF") shows for each case the portion of fractures with r c greater than the given value. In general, the characteristic radii of the matrix blocks formed by the dominant sub-network are larger than those formed by the entire fracture network. The reason is that the dominant sub-network contains fewer of the original fractures [17].

Results for cases with power-law aperture distribution
The other parameter used in this paper to represent the sizes of matrix of varying shapes is the characteristic length which is defined in Eqs. (8) and (9) dominant sub-networks for one realization is presented. Similar results are obtained: for the cases with the broadest power-law aperture distribution, the equivalent length for the dominant sub-network is around twice that on average for the original fracture network. The ratio becomes smaller as the aperture distribution is narrower, and decreases to 1.5 for a = 6.

Results for cases with log-normal aperture distribution
Similar to the results for the cases with power-law aperture distributions, as more fractures are eliminated according to the flowsimulation results, the equivalent matrix-block radius and equivalent matrix-block length become larger for the cases with lognormal aperture distributions (Figs. 7 and 8), and the overall trends are almost unaffected by the breadth of the aperture distribution (Figs. 7d and 8d).
However, the dominant sub-network, defined in our previous study as the sub-network retaining 90% of the original equivalent network permeability after eliminating a portion of fractures, is strongly affected by the aperture distribution: for the broadest aperture distribution cases (r = 0.6), the flow-path length is roughly 20% of the total length of the original fracture network (l b =l o ¼ 0:2), while for the narrowest aperture distribution (r = 0.1), the ratio l b =l o is approximately 0.6 ( Fig. 9). Correspondingly, as presented in Fig. 7, the ratio between the equivalent radius of the dominant sub-network and that of the entire fracture network (r b e =r o e ) is around 3.5 when r = 0.6 (the red dashed line in Fig. 7d), which is the broadest aperture distribution examined here. The ratio r b e =r o e decreases as the aperture distribution becomes narrower, and reaches roughly 1.5 for the cases with the narrowest aperture distribution with r = 0.1 (the red dotted line in Fig. 7d). The equivalent matrix length (L b e ) presents similar behaviour as the equivalent matrix radius (Fig. 8).
Similar to our approach with the cases of power-law aperture distributions, we select one realization for each value of r, and examine the characteristic-radius and characteristic-length distributions of individual matrix blocks in the original fracture network and the dominant sub-network. As presented in Figs. 10 and 11, the characteristic radius and the characteristic length of individual matrix blocks formed by the dominant sub-networks are larger than those of the matrix blocks formed by the entire fracture network.

Discussion
Naturally fractured oil reservoirs, like all reservoirs, are exploited in two stages: primary recovery and secondary recovery. The recovery mechanisms are different in these two processes. During primary production, the reservoir is produced by fluid expansion. The pressure drops rapidly in the fractures because of the high permeability, while, in contrast, the matrix remains at higher pressure. This creates a pressure difference between the fracture and the adjacent matrix block, and in turn, leads to flow of oil from the matrix to the fracture. In this scenario, as long as all the fractures are much more conductive than the matrix, one might expect that all connected fractures are conductive enough to bring oil from the matrix to the wells.
In secondary recovery, since the fractures have much higher permeability than the matrix, water from an injection well invades  the fractures much faster than the matrix. The water rapidly flows through the fracture network and surrounds a matrix block. If the matrix block has water-wet characteristics, water imbibes into the matrix block because of capillary pressure. In many cases, this is a counter-current imbibition process: water and oil flow in opposite directions, although co-current imbibition is faster and can be more efficient than counter-current imbibition [54]. In this case, the cumulative oil recovery from matrix blocks surrounded by water in fractures can be scaled by an exponential equation [55]: where R is the recovery, R 1 is the ultimate cumulative recovery, k is a constant and t is time. Mattax and Kyte [56] redefined the scale equation for imbibition recovery through experimental investigations: where and t D is a dimensionless time, k is the permeability of the matrix block, / is the porosity of the matrix block, r is the interfacial tension, l w is the water viscosity, L m is the matrix block size, and t is imbibition time.
Kazemi et al. [46] further modified the Mattax and Kyte's scaling equation by introducing the shape factor (F s ), which considered the effect of matrix block shapes and boundary conditions. Eq. (13) becomes The definitions of the shape factor (F s ) and the characteristic length (L c ) are described in Eqs. (8) and (9). Thus, when the shape factor of a matrix block becomes smaller, or the characteristic length becomes greater, it takes a longer time to recover a certain portion of oil. As presented in Figs. 6 and 11, in general, for all the cases the characteristic length (L c ) of the matrix blocks formed by the dominant sub-networks are larger than that of the matrix blocks formed by the entire fracture networks. This implies that the rate of oil recovery from the matrix blocks is overestimated if the entire fracture network is considered to take part in the waterflooding process.
Eq. (15) is for one matrix block; here we apply this formula to the average value of a series of matrix blocks in order to estimate the effect of matrix block sizes on the speed of oil recovery during the waterflooding process. For the cases with broadest aperture distributions (a = 1.001 for a power-law aperture distribution and r = 0.6 for a log-normal aperture distribution), the characteristic length (L c ) of the matrix block formed by the dominant sub-networks are on average about twice and three times larger than that formed by the entire fracture networks, respectively (Figs. 5 and 8). This suggests that, for the cases with the broadest power-law aperture distribution, if the entire fracture network is considered to take part in the water-flooding process, the estimated imbibition recovery rate from matrix blocks can be on average four times faster than if only the dominant sub-network is taken into account. For the cases with the broadest log-normal aperture distribution, it can be on average nine times faster. Even for the cases with the narrowest aperture distributions (a = 6 for a power-law aperture distribution and r = 0.1 for a log-normal aperture distribution), the imbibition recovery process can be more than twice as fast.
As discussed above, most water does not flow through the entire fracture network, but mainly through the dominant subnetwork, and the remaining portion of the fractures can be ignored without strongly affecting the overall flow. Considering the recovery mechanism of waterflooding, if the entire fracture network is taken into account, it means we assume the water flows through all the fractures. Then the sizes of matrix blocks formed by fractures and typical fracture spacing are smaller than in the dominant sub-network, which makes the estimated imbibition recovery process faster than it really is. Since dual-porosity/dual-permeability models do not represent fractures explicitly, but assign average properties to grid cells, taking the entire fracture network into account may lead to an inaccurate shape factor, and in turn, give rise to inaccuracy for dual-porosity/dual-permeability simulation of the secondary recovery. The same argument would apply to EOR processes, which depend on penetration of the injected fluid through the fracture network, as waterflooding depends on water.
This suggests that the shape factor for dual-porosity/dualpermeability simulation should depend on the process involved. Specially, it should be different for primary and for secondary or tertiary recovery. For primary recovery, all fractures should be included; for waterflooding or EOR, taking into account only the dominant sub-network which carries almost all water might give a better estimation.

Conclusion
Even in a well-connected fracture network, injected water does not flow through the entire fracture network; it mainly flows through a dominant sub-network which is strongly affected by the aperture distribution. The remaining fractures can be ignored without strongly affecting the overall flow through the fracture network.
The typical fracture spacing and sizes of matrix blocks defined by the entire fracture network are generally larger than those formed by the dominant sub-work which carries most of the flow. If the typical fracture spacing used to calculate the shape factor for a waterflooding process accounts for the entire fracture network, it means the water is assumed to flow through all fractures and all fractures participate in the imbibition process, which is not the case. The shape factor calculated by taking all fractures into account may lead to inaccurate dual-porosity/dual-permeability simulation of the water-flood process. A similar argument applies to EOR; the injected EOR agent does not flow equally through all the fractures.
This suggests that the shape factor for dual-porosity/dualpermeability simulation should be different for primary and for secondary recovery and EOR. Specifically, for primary recovery, all fractures should be included, while for the processes in which delivery of injected fluids plays a limiting role, such as secondary recovery and EOR, the characteristic matrix-block size or shape factor utilised in simulation should be larger.