Impact of Effective Normal Stress on Capillary Pressure in a Single Natural Fracture

Multiphase fluid flow through rock fractures occurs in many reservoir applications such as geological CO 2 storage, Enhanced Geothermal Systems (EGS), nuclear waste disposal, and oil and gas production. However, constitutional relations of capillary pressure versus fluid saturation, particularly considering the change of fracture aperture distributions under various stress conditions, are poorly understood. In this study, we use fracture geometries of naturally-fractured granodiorite cores as input for numerical simulations of two-phase brine displacement by super critical CO 2 under various effective normal stress conditions. The aperture fields are first mapped via photogrammetry, and the effective normal stresses are applied by means of a Fast Fourier Transform (FFT)-based convolution numerical method. Throughout the simulations, the capillary pressure is evaluated from the local aperture. Two approaches to obtain the capillary pressure are used for comparison: either directly using the Young-Laplace equation, or the van Genuchten equation fitted from capillary pressure-saturation relations generated using the pore-occupancy model. Analyses of the resulting CO 2 injection patterns and the breakthrough times enable investigation of the relationships between the effective normal stress, flow channelling and aperture-based capillary pressures. The obtained results assist the evaluation of two-phase flow through fractures in the context of various subsurface applications.


INTRODUCTION
Understanding multiphase fluid flow behaviour in deep fracture-dominated reservoirs is crucial for a variety of applications, such as geological CO2 storage (Ringrose et al. 2013), Enhanced Geothermal Systems (EGS) (Vogler et al. 2017b), nuclear waste disposal (Cassigara et al. 2005), and oil and gas production (Khelifa et al. 2014). When more than one fluid phase is present in the void space of a rock fracture, the different phase will interfere with one another, thus reducing the effective/relative permeability. Additionally, it is important to consider possible discontinuities of some properties at phase boundaries, such as pressure. At the interface between fluid phases, the surface tension effects lead to capillary pressure phenomena, which is defined as the difference between the pressure in the non-wetting phase and the wetting phase. Such capillary pressure and relative permeability functions must be understood in order to properly study multiphase flow in subsurface reservoirs. When laboratory data is not available, it is common to adopt constitutional relations of capillary pressure (Pc) and relative permeabilities (kr) as functions of liquid saturation. Brooks-Corey type and van Genuchten type curves, traditionally developed for porous media, have both been used to describe the capillary behaviour in fractures (Reitsma and Kueper 1994, Yang et al. 2013, Wang et al. 2017. Field-scale numerical simulations often assume zero capillary pressure in fractures, but this might result in underestimation of fluid recovery in deep reservoir systems (de la Porte et al. 2005).
The void space of real rough-walled rock fractures is commonly modelled as a two-dimensional heterogeneous porous medium, where the aperture varies as a function of position in the fracture place and flow occurs in intersecting channels of varying aperture (Pruess and Tsang 1990). This concept is convenient when using traditional porous media numerical simulators, allowing Darcy's law to be applied to fracture flow with only small changes and matrix/fracture interactions to be more easily modelled. In this scenario, some researches build capillary pressure-saturation relations (Pc-S relations) dependent on the aperture distribution, which are suitable for flow in a single fracture. Generating the Pc-S curve includes determining the remaining wetting phase saturation at equilibrium status during drainage, for each condition of applied pressure -which can eventually be associated to the capillary pressure (Wang et al. 2017). In numerical modelling, three criteria are usually considered in analysing the saturation in a mesh grid (hereafter as an element): (i) the capillary pressure criterion, which considers that the non-wetting phase can only invade an element if its pressure overcomes the capillary pressure plus the wetting phase pressure; (ii) the non-wetting phase accessibility criterion, which states that the non-wetting phase can only invade an element if there is a continuous path from upstream to the element; (iii) the wetting phase trapping criterion, which assumes that the wetting phase can only leave an element if there is a continuous path from the current element to a downstream neighbouring element, and otherwise it will be trapped. In consonance with this methodology, numerical percolation models are often applied to calculate the Pc-S relations. One of the earliest models was developed by Broadbent and Hammersley (1957), who used the standard percolation 2 (SP) theory to study capillary flow in fractures. Since then, improvements to this SP model have been proposed, including most importantly the invasion percolation (IP) model, proposed by Wilkinson and Willemsen (1983), in which the trapping of the wetting fluid started to be considered. It is accepted that these models can satisfactorily predict the slow, immiscible displacement processes in narrow fractures (Amundsen et al. 1999).
When modelling two-phase flow through a single fracture as a porous medium, some properties can be determined directly from the local aperture value. In numerical models, the parallel-plate approximation is usually adopted to an element-wise representation of the local fracture aperture (Pruess and Tsang 1990). Here the element area and variation of local apertures inside the element are usually fairly small. It is commonly assumed that the relation of capillary pressure versus fluid saturation must be determine in order to simulate two-phase flow through porous media (Wang et al. 2017). However, for two-phase flow through rough fractures the capillary pressure can be directly calculated using the local aperture value. This implies that the elements are so small that they are either filled with the wetting or the non-wetting phase, and that the Young-Laplace's equation can describe the local capillary pressure. During a drainage process of externally imposed injection, for example, the nonwetting fluid will only invade an element (local parallel-plate opening) if the pressure drop between the two neighbouring elements is greater than the capillary pressure defined by the local aperture. This concept of local aperture-based capillary pressures has been used by various researchers (Kuerper and McWhorter 1991, Reitsma and Keuper 1994, Yang et al. 2013, Huo and Li 2014, and appears more physically sound than the traditional correlations that relate capillary pressure to fluid saturation. Since both relative and absolute fracture permeability depend on the local aperture, flow though fractures further depends on the effective normal confining stress which controls fracture opening and closuring. It is well known that an increasing effective normal stress induces: a decrease of the mean fracture aperture, an increase of the contact area, and a potential increase of tortuosity (Tsang 1984). These changes are tightly coupled to the fracture opening and corresponding transport properties (Pyrak-Nolte et al. 1987). For multiphase flow, both the intrinsic permeability of the fracture and the capillary pressure and relative permeabilities are affected by the effective normal stress, as these properties are also functions of the aperture distribution (Reitsma andKueper 1994, Liu et al. 2013). Huo and Li (2014) have modelled stressdependent capillary pressure curves by using aperture distributions measured in the laboratory, and they concluded that increasing compressive stress results in a decrease in mean aperture and an increase in the variance of the aperture distribution. In contrast, Muralidharan et al. (2004) found that the variance of the aperture distribution decreases with stress, while Liu et al. (2013) observed that it stays constant. Nevertheless, these changes affect the fracture permeability. Therefore, fluid injection/production in fracture-dominated media will inevitably be influenced by variations in effective normal stress (Zhang et al. 2019).
The aperture distributions of rough-walled fractured rock specimens are commonly determined under zerostress conditions, and therefore do not account for stress effects on the fracture aperture distribution (Pyrak-Nolte et al. 1987). Various measurement techniques have been used to characterize aperture fields, such as X-ray microtomography imaging, laser scanning, photogrammetry, and linear profilometer scans (Pruess and Tsang 1990, Caulk et al. 2016. Aside from X-ray tomography scans performed in stressed cores, these techniques give rise to aperture fields under zero effective stress conditions, apparently not applicable to reservoir stress condition. Therefore, in order to obtain representative fracture geometries for studies of hydraulic properties under in situ stress conditions, contact mechanics models should be applied. In recent studies, fast Fourier transform (FFT)based convolution has increasingly been used to solve contact mechanics algorithms, with the advantage of being computationally efficient (Jacobs et al. 2017, Kling et al. 2018. The goals of this study are to (i) map aperture fields of naturally fractured granodiorite specimens via photogrammetry; (ii) numerically obtain aperture fields for different effective normal stresses and build the corresponding Pc-S relations based on the pore occupancy model; (iii) quantify the impact of effective normal stress on two-phase flow of scCO2 and brine via numerical simulations using aperture-based capillary pressures.

METHODS
Naturally fractured granodiorite cores were provided by the Deep Underground Geothermal (DUG) Lab at the Grimsel Test Site (GTS) in Switzerland ( Figure 1) (Tanaka et al. 2014, Vogler et al. 2016. Two laboratory specimens with a diameter of 2.5 cm and length 3.0 cm were chosen for aperture field characterization: N16 (tensile fracture) and N18 (shear fracture). Classification of the fractures followed criteria described by Vogler et al. (2017a).
The fracture surfaces of both specimens have been digitized with a photogrammetric scanner, ATOS Core 3D scanner. The calibration of the photogrammetry scanner was performed in such a way that the length deviation errors were between 0.009 and 0.027 mm, and the optimized calibration deviations were of 0.014 ± 0.001 pixels (Vogler et al. 2016). Further surface alignments allowed determination of the fracture aperture under zero stress condition (Figure 2), as described by Vogler et al. (2017a). Specimen N16, which was a tensile fracture, presented higher mean aperture (185.1 μm) than that (79.4 μm) of Specimen N18, which was a shear fracture.   Once the zero-stress fracture aperture fields were obtained, aperture distributions under non-zero effective normal stresses (2 MPa and 10 MPa) were calculated using a deterministic FFT-based contact mechanics code model developed by Pastewka (Jacobs et al. 2017, Kling et al. 2018. Figure 3 and Figure 4 show the aperture fields of both specimens (N16 and N18), as well as the aperture distributions and the variation of contact area and mean aperture with increasing stress, obtained from the contact model. Higher effective normal stress leads to lower mean aperture of the distribution and more contact area, which yields a less uniform aperture field (with higher coefficient of variation). These different aperture fields were then used to determine the aperture-based capillary pressures and to compose the meshes used for the numerical simulations.

Capillary Pressure Modelling
Two models of capillary pressure were studied: the direct use of Young-Laplace's equation on the numerical simulations and the fitting of van Genuchten (VG) Pc-saturation correlation into curves constructed based on the pore occupancy model.
The difference of the pressures of two phases, on a tensioned surface in equilibrium, can be described by equation [1], known as the Young-Laplace equation: where r1 and r2 are the principal radii of curvature of the interface and γCO2-brine is the surface tension between the CO2 and brine phases (Pruess and Tsang 1990).
When modelling quasi-static displacement within a horizontal, variable aperture fracture, it is generally accepted to use Young-Laplace's equation to calculate the local capillary pressure, which establishes the pressure required by the gaseous phase to occupy an element, under the assumption of a parallel-plate model (Murphy and Thomson 1993). In this scenario, given a local aperture b, equation [1] can then be adapted considering the principal radii of curvature of the interface to be r1 = b/(2cosθ) and r2 = ∞, where θ is the contact angle between the brine meniscus and the fracture wall: A second approach to apply the Young-Laplace's equation in numerical simulations is to compute the capillary pressure-saturation curve (CPSC), based on the pore occupancy model (Liu et al. 2013, Huo andLi 2014). This method relies on the assumption that the local apertures will be filled with either wetting phase or non-wetting phase, assuming small element area and aperture variation within the element (Wang et al. 2017 [3] Considering that a local aperture below bth will be preferentially filled with wetting phase, the saturation of the wetting phase can then be computed as given in equation [4]. Hence, for each value of capillary pressure, a value of water saturation can be computed. where Vw is the total wetting phase volume within the entire fracture, V is the total volume of entire fracture, Ne is the element number within the entire fracture, Δx and Δy are the side lengths of an element, and bi is the aperture of the i-th element. This computation of the CPSC considers only the first criterion for non-wetting phase to invade an element during drainage, based on the standard percolation theory, which is the capillary pressure criterion. The criteria of accessibility and trapping are disregarded during this analysis, what might lead to underestimation of the residual wetting phase saturation during the numerical simulations.
The capillary pressure and the corresponding water saturation were computed at steps of 0.01 kPa, and the results are plotted as capillary pressure-saturation curves, shown in Figure 5. Since a single fracture can be modelled as a two-dimensional porous medium (Wang et al. 2017), the CPSC curves were then fitted using the van Genuchten (VG) correlation (equation where Swr is the residual saturation of the wetting phase. The VG fitting parameters for each stress condition and specimen are listed in Table 2. These parameters were used for the second batch of numerical simulations, using PorousFlow module of MOOSE framework and its built-in VG capillary pressure function, as it is explained in section 2.3.

Numerical Modeling
To compare the impact of effective normal stress on the two-phase flow of brine displacement with scCO2 injection into fractured specimens, a MOOSE framework-based in-house application was used. The fracture was conceptualized as a two-dimensional squared lattice of mesh elements, and the aperture field was taken from the photogrammetry scans of the naturally fractured specimens after the normal stress was numerically applied (Figure 6). Aperture of each element was evaluated using the parallel-plate assumption with a constant opening, and the fracture permeability was calculated using the so-called local cubic law, which is the classical model of permeability for a fracture (Wang et al. 1982, Zimmerman andYeo 2000). The PorousFlow module from MOOSE framework was used, in combination with equations of state for both brine and scCO2, and the two models for capillary pressure mentioned above were considered for each aperture field. The initial conditions were set to 100°C for the temperature and 15 MPa for pressure (see Figure 6). This pressure is kept constant at the outflow boundary. An inlet constant injection flux of 5.6×10 -7 kg/s mimics the conditions of potential future laboratory experiments. Heat transport is not modelled hence the temperature remains constant over time. The left and right edges were considered as no-flow boundaries ( Figure 6). Additionally, the interfacial tension between brine and CO2 was considered constant at 30 mN/m (Pereira et al. 2017) and the contact angle was considered to be zero, i.e., complete wetting.

RESULTS AND DISCUSSIONS
The simulation results show that the effective normal stress has a direct impact on the two-phase flow in single fractures. Under an effective normal stress of 10 MPa, in both specimens, the production of scCO2 occurred earlier than under zero effective normal stress (Figure 7), and more fingering effects were observed (Figure 8). These observations are confirmed by the flow-through laboratory experiments by Raven and Gale (1985), that more tortuous flow channels in natural granite fractures were observed under increasing compressive stress. If these tortuous fingerings are formed in the direction of the flow, which seems to have been the case for the studied specimens, the injection fluid reaches the production side through preferential paths more rapidly. Additionally, the mean cross-sectional area of the fracture decreases with increasing effective normal stress. Hence, with a constant injection mass flow rate throughout all numerical experiments the breakthrough occurs earlier.
The aforementioned observations can be explained by the increased fracture closure caused by higher normal stress, which promotes a decrease in mean value and standard deviation of the aperture distribution (Table 1) and an increase in contact area (Figure 3). For Specimen N18, the increase of 10 MPa in effective stress caused the mean aperture to decrease 50%, and the fractional contact area to increase up to 1.6% (Figure 4). For Specimen N16, even though the compressive stress of 10 MPa yielded a mean aperture only to around 70% of its original value, it was enough to increase the contact area to enhance the channelling of the injection front ( Figure 8). With more contacted asperities under higher stresses, more preferential paths are created and causes an early breakthrough of the injection fluid to reach the outlet boundary and a lower sweep efficiency during the displacement.
A higher normal stress also contributes to higher residual water saturation, as water tends to get trapped in smaller apertures due to capillarity. Figure 9 depicts a comparison between the remaining water saturations (Swr) in Specimen N16 after around 3 hours of injection, during simulations using Young-Laplace's equation and simulations considering capillary pressures equal to zero. The remaining overall water saturation is calculated as the total volume of water remained in the fracture at a certain time, divided by the total volume of the fracture. The remaining saturations increase from 0.16% to almost 0.30% when the stress is increased from 0 to 10 MPa, for the simulations where Young-Laplace's equation was applied. The capillary pressure-saturation curves were also influenced by the effective normal stress ( Figure 5). The influence was more significant in the shear fracture of Specimen N18, in which the initial mean aperture was smaller. Overall, the CPSC became steeper as the effective normal stress increased, suggesting that the capillary behaviour shifted from more fracture-like to more porous media-like. Similar trend in CPSC was also observed by Huo and Li 2014, in numerical simulations using data from flow-through laboratory experiments with increased normal stress. When the aperture field is dominated by higher values, the CPSC is flat, or more fracture-like.
The entry pressure, which is defined as the minimum non-wetting pressure at which the non-wetting phase starts to displace the wetting phase, is associated to minimum aperture values and therefore may also be affected by stress. The decreasing mean aperture caused by increased stress is expected to increase entry pressure. This was the case of Specimen N18, whose entry pressures varied from 0.18 to 0.21 kPa, depending on the applied normal stress ( Figure 5). In contrast, this effect was not observed in Specimen N16, where an entry pressure of around 0.05 kPa was observed for all stress conditions. For Specimen N16, effective normal stresses higher than 10 MPa are necessary to change the entry pressure, due to the higher values of aperture in this tensile fracture.
The simulation results of the models in which van Genuchten (VG) equation was used for Pc and the results of the ones in which the Young-Laplace's equation was applied were similar (Figure 7), suggesting that the VG equation achieved satisfactory modeling of the Pc-S curves of single fracturesat least for the investigated range of flow rate. The Pc curves under different stresses were fitted with different VG parameters (Table 2), which implies that the effect of compressive stress can be modelled by tuning the fitting parameters (a and m) of the VG function. Wang et al. (2017) have correlated the parameter α of VG equation to the reciprocal of the non-wetting phase entry pressure, in the CPSC constructed from their laboratory experiments with a single fracture. In porous media, the parameter α is indeed inversely proportional to the non-wetting fluid entry pressure value, but this conclusion may not necessarily be valid in fractures. In the present work, lower values of a were also associated to Pc-S curves under higher stresses, which are related to higher entry pressures, although a direct reciprocal relationship between a and the entry pressure was not observed. The same trend was observed for parameter m, whose reciprocal can be associated to the slope of the CPSC curve and therefore also to the heterogeneity of the aperture field.
Although the use of van Genuchten porous media capillary pressure function can seem convenient for modelling flow through single fracture, it may not be representative if the aperture fields are significantly heterogeneous, in which case the Pc-S curve has a characteristic of a bimodal curve. The aperture fields studied in the present work were not highly heterogeneous, and therefore the VG equation fit well the curves calculated based on the pore occupancy model (equation [4]). The small "bumps" observed in the curves of Specimen N18 were not large enough to prohibit the modelling using the traditional VG correlation, but this could not have been the case if a higher effective normal stress was applied. Therefore, for the cases in which the fracture aperture fields are sufficiently homogenous, the van Genuchten equation and its fitting parameters can be representative of the CPSC from aperture distribution under different effective normal stresses.

CONCLUSIONS
Two aperture-based capillary pressure models were used to numerically model the injection of supercritical CO2 into brine-saturated single fractures, under different stress conditions. For the simulations, naturally fractured specimens were scanned via photogrammetry and a contact mechanics model was used to numerically obtain aperture fields for different effective normal stresses. The results of the simulations showed that multiphase flow in fractures is significantly affected by the effective normal stress applied to the rock formation. If this stress increases, for example due to the depletion of the subsurface reservoir, the mean value of the aperture distribution and its standard deviation tend to decrease. This directly influences the capillarity of the fractured system, contributing to a more pronounced trapping of wetting phase and fingering during drainage. Additionally, flow through single fractures was accurately represented by both models of capillarity in fractures, which were Young-Laplace's equation and van Genuchten function fitted from curves calculated based on pore occupancy model. This work enables better understanding of applications such as enhanced geothermal systems using CO2 and geological CO2 sequestration in fracture-dominated reservoirs. During long term operation of these systems, the effects of increased effective normal stress on multiphase flow should be a major concern, considering the risk of poor sweep efficiency and/or unintended early breakthrough of the injection fluid.