Does the symmetry of absolute permeability influence relative permeability tensors in naturally fractured rocks?

Relative permeability is a tensorial property directly dependent on absolute permeability tensor which has been taking both symmetric and non-symmetric (full tensor) forms in the literature. Imposing symmetricity during upscaling absolute and effective phase permeabilities could alter both diagonal and off-diagonal terms of relative permeability tensors, specifically in naturally fractured reservoirs (NFRs) where the diagonal components of the permeability tensor could take values with different orders of magnitude. Utilizing a discrete fracture and matrix (DFM) modelling approach, in this paper, we quantify the effects of forcing the symmetricity on absolute permeability tensors on relative permeability tensors of fractured samples in different scales. We also determine the circumstances at which this change in methodology causes a huge difference in the diagonal components of relative permeability tensors.


Introduction
A considerable percentage of the world's conventional oil and gas is reserved in naturally fractured reservoirs (Pirker and Heinemann 2008;Roehl and Choquette 1985;Saidi 1983;Tran et al. 2007). A predictive simulation of these reservoirs requires an accurate fracture and matrix characterization, especially in terms of relative permeability functions which have a significant influence on the two-phase displacement front and sweep efficiency. Thus, a comprehensive understanding of fracture-matrix equivalent relative permeability, k r i , is needed to realize opportunities to enhance recovery from NFRs.
It is difficult to precisely predict water-cut and recovery in naturally fractured reservoirs as flow mechanisms across fractures and matrix are complex (Gilman 2003;Tanaka et al. 2010). Viscous, gravitational, and capillary forces are the main deriving forces determining the ultimate recovery [e.g. Torsaeter et al. (1987) and Sedaghat and Azizmohammadi (2018)]. Capillarity is usually strongly active within the matrix. This leads to strong capillary pressure gradient across the fracture-matrix interface. Viscous and gravitational forces are the dominant forces within fractures. Classical relative permeability models (e.g. Brooks and Corey 1964) assume uniform saturation distribution within the rock. NFRs are highly heterogeneous, and consequently, saturation distribution is not uniform (Behbahani and Blunt 2005;Rangel-German and Kovscek 2006;Uleberg and Kleppe 1996). Hence, classical relative permeability functions are not suitable at all to be applied on NFRs.
Ensemble relative permeability in fractured rocks is strongly dependant on fracture properties, i.e. fracture density, fracture orientation, and fracture connectivity (Fahad et al. 2017). Ensemble relative permeability is also strongly rate dependant. This affects fracture-matrix counter-current imbibition flux by altering the balance between viscous, capillary, and gravitational forces (Matthai et al.

3
2012; Matthai and Nick 2009) which makes the upscaling of multi-phase flows in NFRs very challenging.
Experimental investigation of relative permeability in both small (Bedrikovetsky 2013;Corey et al. 1956;Eichel et al. 2005) and field scale (Rustad et al. 2008;Yeh et al. 1985) shows that relative permeability can be anisotropic due to non-uniform capillary pressure distribution. Numerical analyses show the same too. Pickup and Sorbie (1996) reproduced the results of fine-grid simulations on a coarser grid with phase permeability tensors. Lohne et al. (2006) extended their saturation-pattern-dependent methodology by calculating relative permeability tensors separately in capillary and viscous dominant-layered rock samples for which the layers were aligned with the coordinate system. Later, an analytical example in which relative permeability appears as an asymmetric tensor was derived by Keilegavlen et al. (2012).
Considering relative permeability as a tensor property, it is determined by multiplying the inverse of the absolute permeability matrix to the effective phase permeability matrix (Keilegavlen et al. 2012). Thus, absolute permeability is a very influential property while determining the relative permeability. In fractured rocks, permeability is strongly direction-dependent (Wu et al. 2004), and flow is aligned with the direction in which permeability is higher (Agheshlui et al. 2018b;Azizmohammadi and Matthäi 2017). This anisotropic condition intensifies the effects of absolute permeability on relative permeabilities.
Tensorial analysis of absolute permeability is broadly done in the literature (Durlofsky 1991;Guerillot et al. 1990;He and Durlofsky 2006;Kasap and Lake 1990;Prévost et al. 2005;Wen et al. 2003). Permeability tensor could take two forms: (1) a symmetric tensor in which the off-diagonal terms are identical, k ij = k ji and (2) a nonsymmetric (asymmetric or full) tensor in which the offdiagonal terms are not identical, k ij ≠ k ji (Borisenko 1979;Durlofsky et al. 1994;Gallouët and Guérillot 1991;White and Horne 1987). An asymmetric absolute permeability tensor is theoretically and numerically possible (Ababou 1988;King 1993). The variation formulation of incompressible fluid flow provided by King (1993) requires the work (J) needed to move fluid through a porous medium and P are velocity vector, permeability tensor, viscosity, and fluid pressure, respectively, be minimized. According to Onsager theorem, local permeability, K(x) , is a symmetric tensor. Since ⋅ ⟨∇P⟩ , the variational formulation does not exist based on averaged properties, and consequently, the averaged absolute permeability, ∕⟨∇P⟩ need not be symmetric. However, the symmetricity is usually forced while calculating absolute permeability tensors (Bear 1972;Gelhar and Axness 1983), the rationale of which is given by Onsager (1931). Such a controversy in determining absolute permeability tensor could directly change both diagonal and off-diagonal terms of relative permeability tensor.
To fulfil the symmetric condition while determining equivalent absolute, k E , and effective phase permeability,  In this work, for the first time, we show how forcing the symmetricity on the absolute permeability tensors influences the relative permeability tensors. Here, we aim to quantify the extent of that influence, if there is any, for different fracture geometries in different scales. Here, we do not discuss the superiority of each method. However, we demonstrate the sensitivity of relative permeability tensors on a commonly made assumption while upscaling absolute permeability.
This paper is organized as follows. First, the tensorial analysis of equivalent relative permeability is explained. Second, the outcrop data are presented, and the characterization methods are discussed. Third, the waterflooding simulation scenarios are explained, followed by the results and a discussion.

Methodology
Since this study focuses on naturally fractured reservoirs mapped in outcrop analogues, gravity effect is ignored as the layers are horizontal. Discrete fracture and matrix (DFM) models were built by converting the fracture line tracings into third-order non-uniform rational B-spline curves (NURBS). The resulting geometry was discretized using unstructured grids with lower-dimensional representation of the fractures (Matthai and Nick 2009;Paluszny et al. 2007). We used two-dimensional fracture maps so the built models are two-dimensional. The rock properties and saturation were discretized as piecewise constants on the elements while fluid pressure was discretized as piecewise linear on the nodes and computed with the FEM [e.g. Zienkiewicz et al. (1977)]. The complex systems modelling platform (CSMP++, Matthäi et al. 2007) was used for the computation of pressure and velocity fields.
Two flow problems need to be solved to obtain a full permeability tensor. In the first, we set constant pressure (Dirichlet) boundary conditions on opposite left ( p left ) and right ( p right ) edges of the model, while no-flow (natural) boundary conditions are applied on the others. While the pressure field is solved using the finite-element method, the velocity field is calculated by post-processing the pressure gradients using Darcy's law. Then, the same process is repeated for constant pressures applied to top and bottom boundaries for the second flow problem (Azizmohammadi and Matthäi 2017).
Equation (2) in matrix form can be written as: where superscripts 1 and 2 denote the two flow directions.

Forcing the symmetricity
Solving the system of Eq. 3 gives the components of absolute permeability tensor. This tensor is not necessarily symmetric. To ensure the symmetry of the permeability tensor, the equation k E xy − k E yx = 0 is added to the system (Eq. 3).

Ensemble relative permeability tensor
The fracture-matrix ensemble relative permeability tensor, E ri , for phase i , is defined as where E is ensemble absolute permeability tensor and E e,i is ensemble effective phase permeability tensor for the phase i , which is calculated in a similar manner to that for the E but for the two-phase flow system. Here, the extended Darcy's law for multi-phase flow (Eq. 5 by Muskat and Wyckoff (1937)) is used. So, for each phase, the system of equations (Eq. 3) is solved.
If the symmetricity is imposed on the system for the absolute permeability calculation, it must also be imposed for the effective phase permeability computation. Having calculated the ensemble effective phase permeability, E e,i , the ensemble relative permeability tensor, E ri , for the phase i , is computed from Eq. (4).

Outcrop models
Three natural fracture geometries taken from metre to kilometre scale outcrops are shown in Figs. 1, 2, 3. Table 1 enlists the fracture characteristics of the three outcrop models. The models contain well-interconnected fractures. This implies that the flow dominantly occurs in fractures. Having a wide range of fracture density and length within the models ascertain that a wide spectrum of flow conditions is covered. The models are well-documented in the literature (Belayneh et al. 2006;Matthai et al. 2012), so they were chosen for our research. Fracture aperture distributions were computed based on the far-field stress regime. The subsurface depth is assumed 1.5 km.
To capture fracture geometry, the trace maps of the subhorizontal rock layers were digitized and then discretized using constrained conforming finite-element meshes with lower-dimensional line-element representations of the fractures (Matthäi and Belayneh 2004;Paluszny et al. 2007).

Model configuration
Aperture has been determined for each fracture segment by geomechanical computations (Agheshlui et al. 2018a). Having the aperture determined, local capillary pressure was calculated over the entire model. Table 2 enlists porosity, determined permeability and aperture, connate water and residual oil saturation and relative permeability model for each region of the models. Fluid pressure distributions were computed for every time step with the finite-element method (FEM). Then, having the pressure gradients determined, piecewise constant flow velocities were calculated by postprocessing using the extended Darcy's law (Eq. (5)). These velocities enter the transport equation that is solved with the finite volume method (FVM) in the sequential IMPES, operator-splitting finite-element-centred finite volume (FECFV) scheme (Bazr-Afkan and Matthai 2011; Bazrafkan et al. 2014).
Fracture aperture is derived from geomechanical computations based on Barton et al. (1985) method that considers fracture dilation during shear in response to the far-field stress. Then, the permeability of the fracture segments is computed from the aperture, a , by the parallel plate law (Snow 1969): The capillary entry pressure is calculated from fracture aperture for every fracture segment: In this work, water-oil interfacial tension is assumed constant, = 0.05N∕m , and contact angle, , is set to 0 to establish a strongly water-wet condition. Since our interest is in the bulk flow through the flat-lying reservoir beds in the periphery of wells, we ignore buoyancy forces and inertia effects. This makes the application of linear relative permeability curves to the fractures defendable (cf., Romm (1972); Fourar and Lenormand (2000); Heinemann and Mittermeir (2014)). The Brooks-Corey (1964) relative permeability model was used for the saturation functions assigned to the matrix. Table 2 provides all parameters used in the simulations performed on the outcrop analogue models. Table 3 enlists the waterflooding scenarios applied to each outcrop model. Pressure is fixed at the boundaries meaning that they are open to flow. A constant pressure gradient is maintained across the models during each run (Matthai et al. 2012;Sedaghat et al. 2016a) to mimic unsteady-state core flooding test in the laboratory. All simulations used constant in-and out-flow pressures chosen to emulate pressure gradients that are seen in fractured hydrocarbon reservoirs during production (Table 3) far from the well, where capillary forces are strong. The balance between the forces changes with saturation, the amount of which is reported in Sedaghat et al. (2017). For the sake of simplicity, fluid viscosity and density were assumed to be constant. Water and oil viscosity are 0.001 and 0.0015 Pa s, respectively. The system is treated as isothermal, and the fluids are assumed Newtonian, and incompressible. Given the relatively small size of the models and the moderate variation in pressure within them, this simplification appears defendable.  (Odling et al. 1999)

Results
In this part, first, the impact of forcing the symmetry condition on absolute permeability tensors itself is analysed. Then, its impact on relative permeability tensors is discussed.

Absolute permeability tensors
Tables 4, 5 show the absolute permeability tensor components of the models with and without applying the symmetry condition. For the Kilve model, neither the diagonal terms nor the eigenvalues have changed. However, all permeability tensor components and eigenvalues of the Hornelen model have been influenced by forcing the offdiagonal terms to be equal. Forcing the symmetricity has also made k * xx and k * yy higher and lower, respectively. As a consequence of such a change in the components of the absolute permeability tensor, the minimum and maximum eigenvalues were obtained lower and higher, respectively, representing a completely distinct flow direction and magnitude. For the Arches model, the xx component has not been affected, however, the yy component and min have taken greater values after applying the symmetricity.
In general, forcing symmetricity on the permeability tensor has considerably altered the yy component and min .

3
Also, the off-diagonal terms in the symmetric permeability tensors (Table 5) are closer to the k E xy component of the nonsymmetric permeability tensors (Table 4), specifically for the cases in which k E xx ≫ k E yy . For the models where xx and yy components are generally closer (Hornelen model), k E xx also experiences a noticeable change and the off-diagonal terms of the symmetric tensor are close to the arithmetic average of the off-diagonal terms of the non-symmetric tensor. Figure 4 shows the components of relative permeability tensors for the Kilve model obtained with and without forcing the symmetricity. The words "general" and "forced" stand for not applying and applying the symmetricity, respectively. As it is shown, although forcing symmetricity for the absolute permeability tensor does not affect the diagonal terms and the eigenvalues, k r yx , in significantly influenced for both phases. Here, forcing the    symmetricity causes a decrease in the magnitude of k r yx , although it does not have a noticeable effect on its trend.

Relative permeability tensors
Here, although a more intense counter-current imbibition could be concluded (Sedaghat 2017) if the symmetricity is not applied, the general interpretation of the two-phase flow behaviour does not seem to be hugely affected by the change in the methodology. For the Hornelen model, forcing the symmetricity not only changes all the relative permeability tensor components, it also makes the eigenvalues, the principal values at which the lateral flows get negligible, not match with the diagonal components (Fig. 5). With the increase in saturation, the distance between the curves plotted on the diagonal values (straight to dotted line) increases, and the trend and the magnitude of curves plotted on the off-diagonal values change. Here, in the case of forcing the symmetry condition, in the x-direction, oil flux will be calculated greater than it in the general case. Water flux will take a smaller value if the symmetry is forced. In the y-direction, oil and water fluxes will be calculated smaller and greater comparing to the general case, respectively. For the Arches model, although forcing the symmetricity changes all the relative permeability tensor components, the eigenvalues still match with the diagonal terms (Fig. 6). It also makes the trend of the diagonal terms and the eigenvalues closer to monotonic curves. For the general case, there is a bump for the k r xxi and k r yyi curves in early time which is not usual for relative permeability saturation functions. Different from the previous case, in the x-direction, forcing symmetry leads oil flux to be calculated smaller; however, water flux is calculated higher. In the y-direction, oil and water fluxes are calculated higher and smaller comparing to the general case, respectively. Figure 7 shows the anisotropy ratio (AR) for the models. Applying the symmetry condition does not have a major influence on AR for the Kilve model (Fig. 7a). However, it does have a significant effect on the Hornelen and Arches models. In the Hornelen model, AR is higher for both phases in the case where the symmetricity is applied. For the Arches model, applying the symmetricity causes AR to be lower and higher in low and high saturation values, respectively. The difference made by switching from the conventional relative permeability computation method (Matthai and Nick 2009;Sedaghat et al. 2017) to the tensorial one  with and without applying symmetricity is shown in Fig. 8. For the Kilve model, this methodological difference only affects the AR of the oil phase (Fig. 8a). However, for the Hornelen model, AR is altered for both phases, specifically in high saturation values (Fig. 8b). For Kilve and Hornelen, forcing the symmetricity intensifies the deviation of the tensorial results from the conventional ones. However, this situation becomes totally reversed for the Arches model (Fig. 8c) where forcing the symmetricity makes the extent of that deviation lower and all the curves are decreasing with the saturation, unlike the two other cases.

Discussions
In this paper, constant pressure boundary conditions are used. However, this type of boundary condition is not the best choice. Periodic boundary condition is proved to be a better choice (Durlofsky 1991;Durlofsky 1992;King 1993;Pickup et al. 1992;White and Horne 1987). We did it to just mimic a typical core flooding process. The simulator we used also did not have such capability especially for fractured models where fracture termination points are symmetrically located on the boundary lines. Besides, although running two-dimensional models help us capture the ensemble multi-phase flow behaviour, the simulations had better be performed in 3D in order to come up with a more robust conclusion, especially for naturally fractured rocks where phase segregation is the most dominant flow behaviour in the fractures.
Effective phase permeability is dependent on flow rate (Akin and Demiral 1997;Alizadeh et al. 2007;Nick 2009) andwettability (Batycky et al. 1981;Heaviside et al. 1987;Nguyen et al. 2005;Sedaghat and Azizmohammadi 2019;Sedaghat et al. 2017). In contrast, absolute permeability is independent of those properties. The results of this paper are constrained to the capillary limit as the pressure gradient is too low, and a strongly water-wet condition. Changing the flow rate and the contact angle could alter the results of this paper.
There are also several upscaling methods proposed to determine relative permeability curves. Steady-state upscaling methods provide relative permeability and capillary pressure curves which are only valid over specific flow conditions, i.e. within capillary and viscous limits (Jonoud and Jackson 2008). In this work, we followed an unsteady-state method to both capture the two-phase displacement front behaviour, and mimic the core flooding tests in the laboratory. We also tried to determine the relative permeability at the viscous-capillary equilibrium following the method proposed by Sedaghat et al. (2016b) to see how intense the effects of the symmetricity would be on the relative permeability tensors. However, we could not establish a viscous-capillary equilibrium state within an affordable simulation time. For such complex models where the number of the fractures is high, not all the upscaling methods could be utilized to determine relative Difference between the diagonal components of relative permeability tensor and the one the relative permeability obtained using conventional method permeability and capillary pressure curves with low simulation cost. The representative elementary volume (REV) concept is fundamental for understanding the geological variability of a parameter (Ringrose and Bentley 2015). The concept refers to the scale at which fluctuations in properties approach a constant value (Bear 1972). To have a reliable analysis of flow in porous media, the flow properties must be analysed within its respective representative pore space. In this paper, we showed that the effects of forcing the symmetricity on absolute and relative permeability could be negligible (such as in Kilve), or huge (such as in Hornelen and Arches). Thus, in some cases, adding this extra constraint (symmetricity) could change the size of REV significantly. Assuming symmetric permeability tensors, Azizmohammadi and Matthäi (2017) showed that REV chiefly occurs within 20-40% of the size of the entire model. During their analysis, they had forced the symmetricity on the absolute permeability tensor. Here, we do the same absolute permeability computation over the increasing sample sizes ( l D = sample length model length ) within the model, however, we do not impose the symmetry condition. Figure 9 shows the maximum component of the permeability tensor divided to the matrix permeability, k D = k max k matrix , for different sample sizes. This analysis shows imposing the symmetry condition has no significant impact on REV size. Non-monotonic relative permeability functions have been reported in the literature (Matthäi et al. 2007). This mainly occurs in early time due to having capillary discontinuity at the inflow boundary. However, the reason for having such a case in Fig. 6 is not due to this boundary effect as the methodology is different. In the Arches model, permeability tensor significantly differs even between two neighbour regions (Azizmohammadi and Matthäi 2017). This could dramatically change the multi-phase flow behaviour even within small ranges of saturation. As a result v xx i and v yy i are not always monotonic which leads to non-monotonic relative permeability tensors. It is shown that forcing the symmetricity has made the curves smoother in this case (Fig. 6).
In this paper, we have determined relative permeability tensors with different assumptions for the permeability tensors. However, we have not assessed their reliability to be used for simulation. In order to verify the implemented method and find out if the absolute permeability should be symmetric, all these tensors must be fed into the coarse-grid model to be simulated with the same initial and boundary conditions as the fine-grid model to see if they match. This is the scope of our next paper.

Conclusions
Assuming absolute and effective phase permeability could be either symmetric or asymmetric tensors, equivalent (upscaled) relative permeability tensors were calculated for three DFM models built based on outcrop analogues. The results were compared with the ones obtained from a scalar analysis . The main conclusions of this work are as follows: • In general, forcing symmetricity on the permeability tensor has the most effect on the yy component and also the minimum eigenvalue of the absolute permeability tensor. • If k E xx ≫ k E yy , then k E xy=yx symmetric approaches k E xy asymmetric .
• If the symmetricity is forced, diagonal components of relative permeability tensors will not necessarily match with the eigenvalues of the tensor, regardless of the type of the methodology. • Forcing the symmetricity could lead to overestimation of oil flux in x-direction comparing to the general case. • In Arches model, forcing the symmetricity makes the extent of deviation of relative permeability obtained using a tensorial method with the one obtained using the conventional method lower. • Forcing the symmetricity does not change the REV size.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.