Observations of the Impacts of Millimeter‐ to Centimeter‐Scale Heterogeneities on Relative Permeability and Trapping in Carbonate Rocks

Carbonate rock reservoirs are dominated by heterogeneity across a large and continuous range of spatial scales. We study the impact of heterogeneities on relative permeability and residual trapping for three carbonate rocks selected for their distinct spatial scales of rock texture. The Indiana limestone comprises millimeter‐scale heterogeneities, the Estaillades limestone consists of half‐centimeter‐scale heterogeneities, and the Edwards dolomite includes decimeter‐scale heterogeneity. Along with routine characterization of rock samples, steady‐state N2–deionized water drainage relative permeability measurements are made for each rock at two distinct total flow rates, at least 1 order of magnitude apart. The variation in flow potential across the core results in observations of fluid distribution, core‐average relative permeability, and residual trapping obtained for a range of continuum‐scale capillary number 0.02

continuum-scale capillary number    Δ 0.02 4.2 Δ c c P H N L P . The relative permeability curves for all rocks shift to the right of the water saturation axis with increasing flow potential; the nitrogen relative permeability increases while the water relative permeability decreases. However, the magnitude of the shift depends on the spatial scale of heterogeneity. An inspection of 3D saturation distributions in the cores and estimation of the capillary numbers of flow shows that the rock with the largest heterogeneity is capillary flow dominated throughout the range of injection rates tested; observations in the Indiana and Estaillades carbonates traverse capillary to viscous dominated flow regimes, with commensurate flow rate dependence in the relative permeability. In all cases, residual trapping is poorly described by the Land model.

MANOORKAR ET AL.
It is heterogeneity in multiphase flow characteristics that have a particularly significant impact when upscaled from submeter to field scales . Small-scale heterogeneity in the pore structure results in variation in capillary pressure characteristics which results in saturation heterogeneity, and thus, relative permeability heterogeneity (Chaouche et al., 1994;Egermann & Lenormand, 2005;Huang et al., 1995;Krause & Benson, 2015;Krause et al., 2013;Perrin & Benson, 2010). Pore structure variation can also result in systematic variation in relative permeability characteristics, but it is possible that this heterogeneity is randomly distributed (Zahasky et al., 2020). It is also difficult to characterize in rocks.
A key challenge in modeling flow in heterogeneous rocks is that they lack a representative elementary volume for the multiphase flow properties. Most laboratory observations are obtained using homogeneous rock samples in which measurements of relative permeability are made at a single high flow rate to negate the impact of heterogeneities Pini & Benson, 2013b;Pini & Krevor, 2019). This gives the flow rate-invariant relative permeability, which is also referred to as the viscous limit, intrinsic, or characteristic curve (Krause & Benson, 2015;Rabinovich et al., 2015;Reynolds & Krevor, 2015). However, for rocks where a separation of length scales in heterogeneities cannot be defined, the viscous limit relative permeability is not relevant to field scale modeling. In this case, what is known as an equivalent flow property should be used (Dagon et al., 2013;Jackson et al., 2018). The equivalent relative permeability is a capillary-number-dependent property which must be measured or inferred at the reservoir condition of interest, for example, at the same fluid flow rate with equivalent thermophysical fluid properties. This may be impractical to obtain experimentally.
An approach developed in Jackson et al. (2018), Krause & Benson (2015), and Krause et al. (2013) mitigates this issue by using laboratory observations to create high-resolution numerical models of the rock cores, including parameterization of heterogeneity in the capillary pressure characteristics. These numerical models can then be used to estimate the relevant equivalent flow properties, relative permeability and capillary trapping, needed at the start of an upscaling work flow. While this workflow has been applied to sandstone rocks with layered heterogeneities, it has yet to be demonstrated for the type of multiscale heterogeneities prevalent in carbonate reservoirs (Feng et al., 2019;Hejazi Sayed et al., 2019;Regnet et al., 2019).
In this work, we carry out the experimental stages of this workflow for carbonate rock samples selected for their distinct length scales in rock structures. We evaluate the impact of heterogeneity on relative permeability and residual trapping by making steady-state observations at a range capillary numbers. The flow-rate-dependent equivalent relative permeability and residual trapping characteristics are observed along with the subcore fluid distribution, which can be used to infer heterogeneity in capillary pressure characteristics. With the selected carbonate rock samples, we are able to investigate the impacts of systematic changes to the size scale of the heterogeneities. We are also able to evaluate the impacts of isotropic and nonlayered heterogeneities, distinct from the rock structures typically observed in sandstone. Finally, the observations provided are the information required to subsequently build numerical models that include the upscaled impacts of these heterogeneities.

Rock Samples and Fluids
Three types of cylindrical carbonate rock samples (cores) are used in this study, the Indiana limestone, the Estaillades limestone, and the Edwards dolomite. We selected samples to demonstrate the dependence of flow and trapping on the length scales of heterogeneity in the samples. Indiana limestone is a calcite cemented grainstone made up of fossil fragments and oolites (Churcher et al., 1991) Nitrogen and deionized water are used as the nonwetting phase and wetting phase, respectively. At the experimental pressure (10-12 MPa) and temperature (20°C-21°C), nitrogen is in the supercritical state. Thermophysical properties are estimated from NIST Chemistry Webbook (2018). The density difference between the two fluids, Δρ = ρ w − ρ nw = 887 kg/m 3 , and the interfacial tension γ = 62 mN/m. The viscosity ratio of the fluid pair is μ w /μ nw ≈ 50.

Core-Flooding System
Core-flooding is performed in an elevated pressure and temperature closed loop flow rig as described in Reynolds and Krevor (2015). Dual high-pressure syringe pumps for N 2 and Deionized (DI) water are used to coinject the fluids. Automatic valves for each pair of pumps control the fluid flow and refill the pumps to maintain the continuous circulation of fluids through the core. A two-phase separator is used to separate the phases at the outlet. Fluids are returned to the pumps from the separator via a water return line connected to the base of the separator and N 2 return line to the top of the separator. The rock sample is placed inside an aluminum core holder in a horizontal orientation. A back pressure pump at the outlet side is used to maintain the pore fluid pressure and a confining pump is used to apply an overburden of 4 MPa over the experimental pressure to the core. Pressure is measured at the inlet and outlet face of the core with pressure transducers. All flow lines and pumps are constructed from a corrosion resistant alloy (Hastelloy). Fluid saturations are measured using a medical X-ray CT scanner using voltage of 120 kV, current of 200 mA, and exposure time of 1 s. Scans of 1 mm thickness are taken for the entire length of the core, with a x-y resolution of 512 pixels or 170 µm which gives voxel dimensions of about (0.17 × 0.17 × 1) mm 3 .
All the experiments are carried out with pore fluid temperature of 20°C and pressure of 10 MPa. At the start of each experiment, air is removed from the flow rig by flowing N 2 through all the flow lines at pressures ranging 0.7-1 MPa and released to the atmosphere. A dry scan of the core is taken with a confining pressure of 4 MPa. The pore pressure is then increased with N 2 and a background scan of the core fully saturated with N 2 is obtained at the experimental pressure and temperature. The flow loop is then depressurized to atmospheric pressure and CO 2 is injected to displace N 2 from the core while it is disconnected from the flow loop. Next, deionized water is injected into the loop at the experimental pressure, dissolving CO 2 remaining in the pores of the rock. A scan of the rock sample fully saturated with water is obtained which is used for calculating saturation and porosity. Table 1 summarizes various properties of the samples, including dimensions, average porosity (measured by X-ray CT), and absolute permeability (measured by applying Darcy's law to multirate injection of DI water at experimental conditions). The porosity, ϕ, is calculated with differential images of the core saturated with air and water:

Routine Petrophysical Properties
where CT and I are the X-ray CT scanner attenuation coefficient converted to numerical values in Hounsfield units. The variables CT water and CT dry are values from images of the water-saturated and dry core, respectively. The X-ray instrument is calibrated such that I water = 0 and I air = −1,000. These values can be in the range of 950-1,050 at high pressures inside the core, due to attenuation of the X-ray beam which results in approximately 5% systematic uncertainty in the porosity measurements.
Fluid saturation in the core sample is calculated using two background scans along with experimental scan at each fractional flow rate. The background scans are one scan of the core saturated with N 2 and one scan of the core saturated with water. Nitrogen and water saturation during drainage are calculated using Equation 2: where the subscripts exp, N 2 , and water refer to scan values obtained during coinjection, when the core is saturated with nitrogen, and when the core is saturated with water, respectively. The slice-averaged porosity and saturation are calculated using slice-averaged CT numbers.
The absolute permeability to water is measured after residual trapping is measured in each experiment. The core is flushed with CO 2 to displace trapped N 2 . Permeability is measured by fully saturating the core with water, injecting water at several constant flow rates, and recording the pressure drop across the core. The absolute permeability is inferred from the measurements in the standard way using Darcy's law. The measured permeability decreases in the order Estaillades limestone (140 mD) > Edwards dolomite (46 mD) > Indiana limestone (24 mD) and shows no correlation with porosity. Transverse, slice-averaged porosities of Indiana limestone, Estaillades limestone, and Edwards dolomite samples are shown in Figure 1. All of the cores are heterogeneous.
The permeability varied more than 50% over the course of the experiments for the Indiana limestone sample but remained relatively constant (<10%) for the Estaillades limestone and Edwards dolomite. This change in permeability could suggest fines migration or rock compaction due to the dissolution of carbonate minerals. These uncertainties dominate the uncertainty exhibited in the relative permeability calculations. We have used the absolute permeability measured at the end of the experiment to infer relative permeability.

Relative Permeability and Residual Trapping
Two sets of steady-state relative permeability core-floods are performed on each core. Drainage relative permeability is measured by coinjecting N 2 and deionized water into an initially water-saturated core. The is increased in a stepwise manner from 0 to 1. Fluids are allowed to reach steady state for each fractional flow, that is, when there is a constant pressure differential across the core length and a constant saturation profile. The graphs for pressure drop across the core as a function of time for each fractional flow in Indiana limestone, Estaillades limestone, and Edwards dolomite can be found in supporting information. It is possible that the saturation would continue to evolve over longer time scales than we can observe in the laboratory. We do not observe any variation in saturation distribution over the course of repeated scans obtained during the 30-min interval for a given fractional flow. The multiphase extension to Darcy's law, given by Equation 3, is used to calculate the observed relative permeability of each phase: where k is absolute permeability of the core with length L and cross-sectional area A. The viscosity of fluid is μ at steady-state saturation for a constant pressure drop, ΔP, and continuous injection at flow rate q. The relative permeability, k r,i , of a phase i is taken as a function of saturation of the phase in the core.
Corrections are not made for capillary end effects. While this is standard practice in core analysis, the standard correction approaches are not applicable to rocks with significant heterogeneity due to the combined impacts of core boundaries and rock heterogeneity (Jackson et al., 2018;McPhee et al., 2015). Thus, we report the observed relative permeability which incorporates impacts of both heterogeneity and core sample end effects. An appropriate use of these observations is in the construction of a numerical model of the rock cores, which are able to reproduce the observed behavior (Jackson et al., 2018;Krause & Benson, 2015). This model can then be used to generate appropriate relative permeability characteristics, including the removal of the impacts of boundary effects, for use in flow modeling.
Experiments are performed at high and low total flow rates to evaluate the impact of heterogeneity by varying the balance of viscous and capillary forces (see Section 3.4). Flow rates are selected based on both the permeability and heterogeneity of the cores, as well as the limitations of the experimental setup. For rocks with high absolute permeability such as the Estaillades limestone (140 mD), 20 mL min −1 is used for the high flow rate observations. In contrast, for the low permeability Indiana limestone (24 mD) and Edwards dolomite (46 mD), 5 mL min −1 is used. The lower flow rate experiments are carried out at 0.5 mL min −1 . This range of flow rates is constrained on the low end by the smallest fractional flow of the nonwetting phase within the pump specifications and on the high end by the pressure differential across the core relative to the confining pressure.
Residual trapping is characterized as the final stage in the same core-flood test, after (  2 1 N f ). We follow the approach of Niu et al. (2015). The capillary pressure gradient at  2 1 N f results in a saturation gradient along the length of the core. This can be used in residual trapping experiments to rapidly create an initial-residual curve. Imbibition is then performed with deionized water at the same flow rate as the drainage component of the experiment (shown in Table 1), until the residual saturation is reached. Thus, a single core-flood is used to obtain a region of the initial-residual curve.
Residual trapping is represented by plotting the initial-residual saturation curve. The Land trapping coefficient, C, is found by fitting Equation 4 to a plot of the residual saturation, , 2 S N r versus the initial saturation, , 2 S N i (Land, 1968

Capillary Pressure Characteristic Heterogeneity
The three-dimensional imagery of saturation in the core is the basis for the characterization of capillary pressure heterogeneity (Egermann & Lenormand, 2005). The role of capillary heterogeneity to permeability heterogeneity is characterized through a dimensionless capillary number describing the ratio of viscous driven flow to the flow driven by gradients in capillary pressure (Virnovsky et al., 2004;Yokoyama & Lake, 1981;D. Zhou et al., 1997). In our experiments, we cover a range of capillary number by changing the flow rates. As the fractional flow rate of N 2 increases in the drainage cycle, the differential pressure across the length of the core decreases because of low viscosity of N 2 relative to DI water. We used the number as defined by Virnovsky et al. (2004), applicable to continuum properties. The capillary number is defined as where H is a characteristic distance between heterogeneous layers and value is given in Table 1 (Hejazi Sayed et al., 2019). It is obtained by inspecting the three-dimensional map of the dry core. The characteristic difference in capillary pressure between layers is given as ΔP c . Different approaches are used for different rock samples to estimate the value. For the Indiana limestone, the pressure differential across the rock core at which low permeability layers are invaded by N 2 is assumed to bound ΔP c . For the Estaillades limestone, there are no clear layers and ΔP c is approximated through the distribution in saturation in a given slice . For the Edwards dolomite, N 2 never invades parts of the impermeable vuggy structure on the downstream section of the core. Hence, the maximum pressure differential observed during the core-flood is taken as ΔP c , although the true value could be higher. These values are tabulated in Table 1.
This is to say that estimates of N c are necessarily subjective and can only really be used to show qualitative trends with varying N c within a rock structure (Virnovsky et al., 2004). Some studies have shown that the transition between capillary and viscous dominated flow regimes might take place in a range of capillary number 1 < N c < 100 for sandstones (Krause & Benson, 2015;Reynolds & Krevor, 2015;Virnovsky et al., 2004). The transition to the viscous limit will be different for different rocks and depend additionally on the geometry of heterogeneity in rocks relative to the flow direction.

Indiana Limestone: Smallest Heterogeneity
The Indiana limestone carbonate rock has the smallest scale of heterogeneity of the rock samples tested, around a millimeter ( Figure 1 and Table 1). The results of the drainage relative permeability tests are shown in Figure 2. At high flow rates, the relative permeability to N 2 ,   for the saturation range of 0.93 ≤ S w ≤ 0.77. The crossover point, where the same value of relative permeability to N 2 and water is achieved, occurs at higher water saturation (S w = 0.78) than at the high flow rate (S w = 0.54).
The increased impact of heterogeneity at the low flow rate serves to increase the mobility of the nonwetting phase and to decrease the mobility of the wetting phase ( Figure 2). The variation in the relative permeability reflects a transition from a state in which heterogeneities control the fluid distribution at low flow potential, to a state in which they are relatively unimportant at high flow potential. It appears that the high flow rate observations approach the viscous limit regime.
The impact of capillary number on the fluid distribution is shown in Figure 3. At the beginning of the drainage displacement, the capillary number is high due to the high water fractional flow. The fluid saturation distribution is homogeneous. As the N 2 fractional flow increases, the capillary number decreases, mainly because of a lower value of pressure differential across the core, associated with the lower viscosity of N 2 . At the lowest capillary numbers, in the low flow rate experiments, layering is observed in all of the images. This shows that a significant transition in the flow regime occurs, from a viscous dominated regime at the initiation of drainage, to capillary dominated at low capillary number. The sectional maps of water MANOORKAR ET AL.
10.1029/2020WR028597 6 of 15 saturation at the center slice along the direction of the flow are shown in Figure 4 (the saturation maps for additional fractional flow are provided in supporting information). At low flow rate, 0.5 mL min −1 experiments, it is evident that the nitrogen is trapped in the layers parallel to the flow direction where as these layers are disappeared for high flow rate, 5 mL min −1 experiments because the entry pressure barrier for new pores is crossed. The threshold capillary number (N c ) for the transition from capillary dominated flow to viscous dominated flow lies between 0.1 and 0.27 which depends on the type of the rock and orientation of the heterogeneity .
A wide range of residual trapping is observed for the Indiana limestone rock sample (Figure 2 (right)). However, the Land model is a poor fit for the high flow rate observations. There is not a systematic impact of the flow rate on trapping.

Estaillades Limestone: Middle Heterogeneity
The Estaillades limestone carbonate core comprises heterogeneities of approximately one half centimeter length scales (Figure 1 and Table 1). The results of the drainage relative permeability tests are shown in Figure 5. At high flow rates, the relative permeability to N 2 ,   for the saturation range of 0.95 ≤ S w ≤ 0.78. The crossover point, where the same value of relative permeability to N 2 and water is achieved, occurs at higher water saturation (S w = 0.78) than at the high flow rate (S w = 0.7).
As with the Indiana limestone, the relative permeability shifts rightward at low capillary number, with higher mobility for the nonwetting phase and lower mobility for the wetting phase ( Figure 5). In this case, the magnitude of the transition between viscous and capillary dominated flow systems is less than for the Indiana limestone core. This manifests as less impact of the flow rate variation on the observed relative MANOORKAR ET AL.   permeability. At the highest capillary number, it appears that the experimental conditions are still farther from the viscous limit.
Residual trapping observations show a higher average trapping of N 2 for low flow rates. It is possible that this reflects a larger scale of capillary trapping than that associated with ganglia snapoff in pores. This mode of trapping involves fluid immobilized upstream of capillary barriers (Krevor et al., 2011). The larger difference in residual trapping for low flow rate and high flow rate experiments for Estaillades limestone is also due to the higher imbibition flow rate as compared to other carbonate cores in this study.
The impact of heterogeneity on the fluid saturation distribution is shown for a range of capillary numbers, N c , in Figure 6. The patchy distribution of fluids does not significantly change across the range of observed capillary numbers suggesting less of a transition between capillary and viscous dominated flow than for the Indiana. The sectional saturation maps at the center slice along the direction of the flow for different fractional flows are shown in Figure 7 which do not show qualitative change in the fluid distribution as compared to Indiana limestone (the saturation maps for additional fractional flow are provided in supporting information).

Edwards Dolomite: Largest Heterogeneity
The Edwards dolomite sample has multiple scales of heterogeneity, with a large impermeable structure near the outlet of the core spanning a few centimeters. Interestingly, the shape of relative permeability curves is similar for both the high and low flow rate experiments with a slight increase in the nonwetting phase mobility at low capillary number (Figure 8 (left)). This is unlike the observations in the other carbonate samples which showed clear flow rate dependencies. As will be discussed below, in this case, the relatively similar behavior at the varying flow rates is because of the relative strength of the heterogeneity. The flow behavior is strongly controlled by the heterogeneity across all of the flow rates.
The residual trapping observations exhibit consistent behavior as the relative permeability. The trapping is on average greater at lower flow rate but is less strongly affected by flow rate than for the other rock samples (Figure 8  The impact of heterogeneity on the fluid distribution is shown in Figure 9. The fluid distribution is highly heterogeneous for both high and low flow rate observations, indicating that the heterogeneity dominates the flow for a wide range of capillary numbers. The sectional saturation maps at the center slice along the direction of the flow for different fractional flows are shown in Figure 10 (the saturation maps for additional fractional flow are provided in supporting information). The fluid distribution does not change qualitatively for low flow rate and high flow rate set of experiments. This is why the relative permeability and trapping are not as strongly effected as with the other rocks.

Conclusions
In this work, we have shown the interplay of different scales of heterogeneity, and capillary number, on observations of subcore fluid distribution, relative permeability, and residual trapping. For all three rocks, the relative permeability curves shift to the right with decreasing capillary number. This is an increase in MANOORKAR ET AL.   the mobility of the nonwetting phase and a decrease in the mobility of the wetting phase. This reflects an increased control of the distribution of fluids by capillary forces which allows the nonwetting phase to flow through high saturation channels (Reynolds, 2016). The extent of the shift is strongly dependent on the strength and length scales of heterogeneity in the rocks. There is the greatest transition, from viscous to capillary dominated flow, in the Indiana limestone, the rock with the smallest heterogeneity. There is very little shift in the relative permeability in the Edwards, where the large heterogeneity imposes capillary controls on the fluid distribution across the range of flow rates observed.
The residual trapping behavior is more complicated. There is a wide range of residual trapping behavior for three carbonate rock types in this study. There could be additional trapping upstream of local capillary barriers for the low flow rate experiments. However, there is significant scatter in the trapping data, and the Land model is not generally a good fit.
MANOORKAR ET AL.   Fluid saturation observed within the cores are consistent with transitions between viscous dominated and capillary dominated flow regimes. In the Indiana carbonate, the observations where the largest transition is observed, there is a clear layering of fluid flow at low capillary number giving rise to a more homogeneous distribution of fluids at high capillary number. In contrast, the fluid distribution remains relatively static across flow rates for the Edwards dolomite.
The experimental observations presented in this paper can be used to construct a core-scale digital rock model to derive equivalent flow properties from the heterogeneous rocks (Hejazi Sayed et al., 2019;Jackson et al., 2018;Krause et al., 2013). This is the first step in an upscaling workflow for reservoir characterization that incorporates the impacts of small-scale heterogeneities in multiphase flow properties.

Appendix: Tabular Relative Permeability Data
Tables A1-A3 are tabular data required to reproduce the graphs in Figures 2, 5,