Characterizing Drainage Multiphase Flow in Heterogeneous Sandstones

In this work, we analyze the characterization of drainage multiphase flow properties on heterogeneous rock cores using a rich experimental data set and mm‐m scale numerical simulations. Along with routine multiphase flow properties, 3‐D submeter scale capillary pressure heterogeneity is characterized by combining experimental observations and numerical calibration, resulting in a 3‐D numerical model of the rock core. The uniqueness and predictive capability of the numerical models are evaluated by accurately predicting the experimentally measured relative permeability of N2—DI water and CO2—brine systems in two distinct sandstone rock cores across multiple fractional flow regimes and total flow rates. The numerical models are used to derive equivalent relative permeabilities, which are upscaled functions incorporating the effects of submeter scale capillary pressure. The functions are obtained across capillary numbers which span four orders of magnitude, representative of the range of flow regimes that occur in subsurface CO2 injection. Removal of experimental boundary artifacts allows the derivation of equivalent functions which are characteristic of the continuous subsurface. We also demonstrate how heterogeneities can be reorientated and restructured to efficiently estimate flow properties in rock orientations differing from the original core sample. This analysis shows how combined experimental and numerical characterization of rock samples can be used to derive equivalent flow properties from heterogeneous rocks.


Introduction
Fluid flow in the subsurface is dominated by rock heterogeneity with length scales ranging from nanometers to kilometers (Ringrose & Bentley, 2015). Heterogeneities at the submeter scale, such as laminae (mm) and bedding (cm), have a large impact on single and multiphase fluid flow properties, and are considered to be an important control on flow at larger scales (Corbett et al., 1992;Li & Benson, 2015;Nordahl et al., 2005;Ringrose & Corbett, 1994). Incorporating accurate heterogeneous flow properties into reservoir characterization is therefore crucial to the successful modeling and prediction of fluid flow in the subsurface , which is important for applications such as groundwater hydrology, petroleum engineering, and carbon sequestration (Ringrose & Bentley, 2015).
A key component in the characterization effort is the sampling of rock cores directly from the subsurface and derivation of flow properties at the mm-m scale in the laboratory (McPhee et al., 2015). For single phase flow, core analysis typically involves the evaluation of porosity, absolute permeability, resistivity, and compressibility. When multiple fluid phases are involved, capillary pressure, relative permeability, and resistivity index relationships are also characterized (Bachu & Bennion, 2008;Dumore & Schols, 1974;Honarpour et al., 1986;Krevor et al., 2012;Pini & Benson, 2013a). Since the multiphase flow properties depend on the saturation and the direction of change in saturation, the hysteretic nature of the functions, and the trapping of fluid phases are also evaluated (Akbarabadi & Piri, 2012;Fulcher et al., 1985;Oak, 1980;Pentland et al., 2011;Pini & Benson, 2017).
The theoretical justification for the use of flow properties derived from rock core analysis is based in Effective Medium Theory; see Renard and de Marsily (1997) and more recently Ringrose and Bentley (2015) for reviews. Measurements of flow properties must be made in a representative elementary volume (REV) of material, whereby the length scale over which the macroscopic properties are derived is larger than the length scale of any heterogeneity within the REV. We draw the distinction here between effective properties, which are generally upscaled quantities where the heterogeneity is much smaller than the averaging scale and equivalent properties, where this separation of scales cannot be strictly defined (Dagan et al., 2013;Durlofsky, 1991). Equivalent properties are both flow rate and heterogeneity dependent, and are more accurate in modeling large scale flow in the reservoir (Renard & de Marsily, 1997).
For multiphase flow properties, the REV for the system is not always clear (Armstrong et al., 2014), especially in the presence of laminae (mm) or bedding (cm) scale heterogeneities (Corbett et al., 1992). In the presence of large heterogeneities, i.e., approaching the scale of observation, effective multiphase flow properties cannot be obtained (Renard & de Marsily, 1997). As a result, most laboratory observations of multiphase flow properties are required to be obtained from homogeneous and isotropic rock samples. Measurements of steady state relative permeability are often performed at a single, high flow rate (typically 10-40 mL min 21 , Krevor et al., 2012;Pini & Benson, 2013a) so that homogeneous saturations are obtained throughout the core. This provides a measurement of the rate-invariant relative permeability in the viscous limit.
The viscous limit relative permeability is often referred to as ''intrinsic'' (Reynolds & Krevor, 2015), ''characteristic'' (Krause & Benson, 2015), or as a ''rock curve'' (Pickup & Stephen, 2000), due to its association with universality. However, as we will demonstrate in this work, the characterization of a rock core solely with this viscous limit relative permeability has little relevance for the description of subsurface flow in heterogeneous rocks. It should be replaced with equivalent relative permeabilities derived on representative rocks at appropriate reservoir conditions.
Deriving equivalent functions experimentally, although possible, places severe demands on the experimental observation. The experiment must be performed with flow rates, fluid pairs, fluid properties, and the orientation of rock heterogeneity with respect to flow as they are in the reservoir (Reynolds & Krevor, 2015). It is not practical to derive these functions across a large range of conditions; as a result no attempts have been made to develop a protocol for effectively measuring properties on heterogeneous rocks. In this work, we show that advances in both core flood experiments and numerical simulation can be combined to overcome the problems associated with deriving equivalent properties on heterogeneous rocks.
A number of studies have shown how saturation heterogeneity, and thus relative permeability heterogeneity, will arise in rocks with small variations in capillary pressure characteristics (Egermann & Lenormand, 2005;Honarpour et al., 1995;Huang et al., 1995;Krause & Benson, 2015;Ringrose et al., 1993). These heterogeneities are often small enough (i.e., DP e < 1 kPa) for the core to be considered homogeneous when evaluated with single phase tracer tests (Reynolds & Krevor, 2015). Although initially focused on oil production, the study of capillary pressure heterogeneity has more recently shifted to understanding gas-water systems associated with CO 2 storage. In these systems, the impact of capillary pressure heterogeneity has been demonstrated experimentally at both the core scale, i.e., Perrin and Benson (2010), Krevor et al. (2011), Shi et al. (2011), Wei et al. (2014) and at intermediate scales below seismic resolution (1-10 m) using meter-scale experiments (Trevisan et al., 2017b). Here it was explicitly shown how layered capillary heterogeneities could increase nonwetting phase trapping by up to 15% in some cases and lead to plume migration time scales of the order of months as opposed to weeks in the homogeneous case.
Numerical simulations at the cm-m scale have demonstrated the larger scale impacts of local capillary pressure heterogeneity, which can immobilize a buoyant CO 2 plume under limited CO 2 supply (Li & Benson, 2015), channel CO 2 through a small fraction of the domain (Meckel et al., 2015) and significantly increase the overall trapping of CO 2 compared to residual trapping alone (Saadatpoor et al., 2010). Furthermore, at the reservoir scale, only capillary flow models can effectively describe the Sleipner field CO 2 plume evolution in the northern North Sea (Cavanagh & Nazarian, 2014), and early breakthrough at the Frio injection site is thought to be largely due to preferential capillary pressure pathways (Hovorka et al., 2006). It is therefore clear that capillary pressure heterogeneities, occurring at length scales of mm-cm, significantly impact the flow and trapping of CO 2 from the laboratory to the field scale, and need to be accurately characterized and incorporated into subsurface modeling efforts.
Incorporating the impacts of small-scale capillary heterogeneity into reservoir scale simulations requires upscaling to be performed, due to the limitations of grid size in numerical modeling. When capillary pressure heterogeneity is incorporated into a dynamic upscaling procedure (as opposed to using simpler capillary or viscous limit upscaling), the resulting equivalent relative permeability is both flow rate and heterogeneity dependent, i.e., dependent on the capillary number (Odsaeter et al., 2015;Rabinovich et al., 2015;Virnovsky et al., 2004). Making use of such a modeling framework, however, requires the development of laboratory based characterization protocols from which the properties of capillary heterogeneity, and their impacts on multiphase flow are initially derived.
Recent work has demonstrated that fluid saturation heterogeneity observed in rock cores using X-ray CT imaging can be used as a sensitive measure to characterize capillary pressure heterogeneity at the REV scale (Egermann & Lenormand, 2005;Krause et al., 2013;Pini & Benson, 2013b;Pini et al., 2012;Reynolds & Krevor, 2015). We build on these experimental approaches and combine this with recently developed modeling techniques, primarily those of Krause (2012), Krause et al. (2013), Krause and Benson (2015), to create a numerical representation of the rock core. The numerical models are then used to derive equivalent relative permeabilities at a range of conditions that are relevant to the subsurface.
Under this paradigm, the goal of the characterization shifts. Conventionally, the focus of core analysis is on obtaining the most accurate viscous limit relative permeability curve directly from experimental core floods (i.e., at very high flow rates in homogeneous cores). Recent modeling efforts have similarly focused on deriving the viscous limit relative permeability from observations on heterogeneous rocks (Krause & Benson, 2015). In the approach presented here, the goal becomes to not only characterize viscous limit properties, but also to derive equivalent properties at relevant reservoir and flow conditions, which are suitable for upscaling and field-scale reservoir simulation.
The characterization presented herein focuses on drainage multiphase flow, which is key to the successful modeling of primary plume migration during CO 2 injection in subsurface aquifers, and in primary hydrocarbon recovery (oil/gas displacing water). To model imbibition processes, typically encountered in secondary hydrocarbon recovery (water displacing oil/gas) or at the trailing edge of a CO 2 plume, incorporating hysteresis and trapping into the flow functions is critical. Recent work by Pini and Benson (2017) has shown how imbibition capillary pressure heterogeneity can be effectively derived, and several models exists for characterizing imbibition relative permeability based on trapping data and/or bounding curves (Killough, 1976;Land, 1968). However, characterizing imbibition processes at varying capillary number in heterogeneous systems introduces significant extra complexity in determining the end point saturations and the path dependence of each phase, which is beyond the scope of the current work.
In the following, we characterize the drainage multiphase flow properties of Bentheimer and Bunter sandstone cores with differing principal orientations of capillary heterogeneity. Routine properties such as porosity, absolute permeability, and viscous limit relative permeability are first derived in section 2.1. Capillary pressure heterogeneity is then characterized in 3-D through the calibration of a numerical model of the rock core with X-ray imagery of fluid saturation data obtained during steady state core flooding experiments. The numerical models are used to predict the equivalent relative permeabilities across a range of capillary numbers, and compared directly to experimental measurements in section 3. Finally, we use the numerical models to predict equivalent relative permeabilities without the constraints imposed by the experimental observations: we remove boundary artifacts, control the structure/orientation of the heterogeneity and perform numerical experiments under conditions impractical to produce in the laboratory.

Characterization of Multiphase Flow Properties on Heterogeneous Cores
Two sandstone cores are used to demonstrate the characterization approach, with properties detailed in Table 1. We use a Bentheimer sandstone with a simple heterogeneity parallel to the axis of flow as described in Reynolds and Krevor (2015), and a Bunter sandstone with primary heterogeneities perpendicular to the flow as described in (Reynolds et al., 2018). The Bentheimer sandstone is a shallow marine sandstone which forms the unit for oil reservoirs in the Netherlands and Germany. It is frequently used in petrophysical studies due to its availability in quarries, high permeability, and homogeneity (Peksa et al., 2015). The sample was composed of 95% fine to medium-grained quartz with a well-sorted grain size distribution and minor feldspar and clays (5%).
The Bunter sandstone sample is a medium-grained sandstone composed mainly of subangular to subrounded quartz grains with a minor component of detrital K-feldspar, clay, and carbonate clasts with a wellsorted grain size distribution. The core sample was obtained from the British Geological Survey Core Store, drilled from the Cleethorpes-1 geothermal borehole at 1,312 m depth.
Since this work is motivated by the impact of submeter scale rock heterogeneity on multiphase flow characterization, we focus on gas-liquid systems at elevated pressures and temperatures, equivalent to those in typical subsurface reservoirs suitable for carbon sequestration . In these systems, low flow potential gas plumes are significantly impacted by capillary pressure heterogeneity (Saadatpoor et al., 2010), which can result in equivalent relative permeabilities which are far from the viscous limit values (Li & Benson, 2015;Rabinovich et al., 2015). We use both Nitrogen-DI water and CO 2 -Brine fluid pairs, so that the impact of capillary pressure heterogeneity can be clearly demonstrated independently of variations caused by the choice of fluid pair. In the DI water experiments with Bentheimer, clay mobilization was monitored through the pressure drop. Repeat experiments on the same core yielded similar absolute permeabilities, indicating no significant swelling (Reynolds & Krevor, 2015).
The importance of capillary heterogeneity in controlling fluid flow depends on the ratio of the viscous or bouyant flow potential to gradients in capillary pressure induced by rock heterogeneity, which can be described with a macroscopic capillary number. In this work, we use the macroscopic capillary number N c presented in Virnovsky et al. (2004) to guide our characterization efforts, where DP is the pressure drop over a macroscopic length scale L and DP c is a characteristic difference in capillary pressure between regions separated by heterogeneity length H. For the systems considered in this work, DP c can be estimated using the difference in entry pressure between the high and low entry pressure regions, with corresponding heterogeneity length scale H. With this and the measured pressure drops, the capillary number can be estimated for the core floods.
This capillary number indicates the transition from a capillary dominated regime (low N c ), where the fluid saturation is controlled by capillary pressure heterogeneity, to a viscous dominated regime (high N c ) where the fluid saturation is primarily controlled by the distribution of absolute and relative permeabilities. The more widely used capillary number n c 5vl=c is indicative of the average balance of viscous and capillary forces at the scale of individual pores. In the absence of chemical surfactants, even corefloods at flow rates much higher than those found in the subsurface exhibit n c < 10 26 , and capillarity is dominant at the scale of individual pores. It is the core scale capillary number N c that provides useful information about the impact of rock heterogeneity.
To characterize the viscous limit relative and absolute permeability, we require N c to be large, in a regime nearing the viscous limit. This ensures that the core has a near homogeneous, 1-D saturation profile through the length of the core, and the relative permeability can be accurately derived without prior knowledge of the 3-D capillary pressure heterogeneity. In the high flow rate regimes considered here for the Bunter and Bentheimer, the maximum N c across the fractional flow regime is of the order 10 1 210 2 , respectively (see section 3.2), which is similar to the viscous limit transition region identified by Reynolds and Krevor (2015) and Reynolds et al. (2018) of 10 < N c < 75. In this viscous dominated regime, the slice-average saturations are largely homogeneous; see Figure S1 in the supporting information for 1-D slice-average saturation maps.
To characterize the capillary pressure heterogeneity, we use observations obtained at flow regimes far from the viscous limit, where the capillary pressure is the main control on the saturation distribution (Pini & Benson, 2017). For the low flow rate Bunter and Bentheimer experiments, the maximum N c across the fractional flow regime is of the order 10 21 210 0 , below the capillary-viscous transition region identified in Reynolds and Krevor (2015) and Reynolds et al. (2018).
The specific flow rates for the steady state core flood experiments are chosen based on the above N c criteria with some prior testing on the core. In practice, for these highly permeable cores (%2 Darcy), we use Q tot in the range 20-40 mL min 21 for the high flow rate experiments, and Q tot in the range 0.2-7 mL min 21 for the low flow rate experiments.
In the following sections, we present the combined experimental and numerical characterization methods, moving from routine rock properties through to capillary pressure heterogeneity and numerical calibration.

Routine Characterization
Here we describe the characterization of routine properties; porosity, absolute permeability, and viscous limit relative permeability. Steady state core flood experiments are performed using the experimental technique described in detail in Reynolds (2016), Reynolds and Krevor (2015), and Reynolds et al. (2018). The core-average experimental data and 1-D saturation profiles used throughout this work are provided in the supporting information. At each fractional flow in the experiments, steady state was obtained once the pressure drop DP and average saturation S w in the core had stabilized (typically DP=Dt < 2kPa/h and DS w = Dt < 0.01/h), see the supporting information and Reynolds (2016) for further discussion.
3-D porosity maps are calculated using medical X-ray CT scanning of the brine saturated rock core through differencing of the saturated image with a dry-scan image, in the standard method described in Akin and Kovscek (2003). For the 3-D porosity and saturation maps, scans were repeated 3 times at each step to increase precision; see section 3 for a discussion on the error quantification for voxel, slice and coreaveraged porosity/saturation. Transverse, slice-averaged porosity maps of the Bentheimer and Bunter cores are shown in Figure 1a, highlighting the homogeneous axial structure of the Bentheimer, and the heterogeneous layering present in the Bunter.
The absolute permeability to brine is calculated at multiple flow rates using the 1-D form of Darcy's law applied over the length of the rock core, shown in Table 1. In this characterization effort, we do not seek to find the 3-D permeability map as in previous studies (Krause & Benson, 2015;Krause et al., 2013), since at the submeter scale capillary pressure characteristics are the primary control on the saturation distribution (Kuo & Benson, 2013).
The high flow rate steady state core flood experiments are used to characterize the viscous limit relative permeability; see Table 1 for a summary of the experimental conditions. To determine a suitable functional form of relative permeability, the experimental slice-average saturation data and pressure drops at each incremental gas fractional flow are used in a numerical regression scheme in the 1-D simulator SENDRA TM . Darcy's law combined with conservation of mass is solved with a given relative permeability parameter estimation. With each iterative solution of the system of equations, the relative permeability function, k rj ðS j Þ is varied using a least squares numerical regression to best fit the experimental 1-D saturation profiles and pressure drops with the simulated values. A P c ðS w Þ function determined from mercury intrusion capillary pressure (MICP) testing (described in section 2.2) is also included in the simulation to take into account capillary end effects. Since the fluid saturations in the high flow rate experiments are effectively homogeneous 1-D profiles in the central region of the core (see supporting information Figure S1), the 1-D numerical regression scheme incorporating capillary pressure can accurately derive the viscous limit (''intrinsic'') relative permeability.
In this work, we find the Chierici relative permeability functional form the most suitable fit for our experimental data (Chierici, 1984). The function is defined by the following equations: where subscript r; g; w; c, and irr refer to relative, gas, water, critical, and irreducible, respectively. The Chierici formulation allows regression of the parameters A, B, M, and L to control the shape and curvature of the relative permeability function. Based on the experimental observations of the saturation maps, we set S gc 5 0 for both cores, and S wirr 5 0.0812 for the Bentheimer and S wirr 5 0.0821 for the Bunter, discussed in more detail in section 2.2. We also set k rg ðS wirr Þ5k rw ðS gc Þ51 for both cores. In this work, we find that k rg ðS wirr Þ has to be high to model the steep increase in k rg as S w is lowered (see Figure 1b), and hence we set k rg ðS wirr Þ51 to obtain the best match with experimental data.
Upon determination of the best fit relative permeability parameters, some manual adjustments are also performed to ensure the curves near the endpoints (which are not reached in the slice average, viscous limit) accurately reproduce the voxel saturations in 3-D simulations at lower flow rates. This generally only involves adjusting the gas relative permeability to allow high water saturations to be reached, i.e., increasing B in equation (2). The viscous limit relative permeabilities for each core calculated through numerical regression are shown in Figure 1b with the Chierici regression fit parameters detailed in Table 1. The discrete experimental relative permeability k rj of each phase j in Figure 1b are calculated using core-average saturation data and pressure drops, i.e.: where q j is the volumetric injection rate, L is the core length, and A is the core area. Note, we assume here that in the viscous limit, the phase permeabilities are equal at the inlet and outlet end caps, i.e., P c 5 0. Equation (3) generates an equivalent relative permeability, which is an upscaled function incorporating the impact of subcore scale capillary pressure heterogeneity and any impact of capillary end effects at the inlet and outlet boundaries of the rock sample. It is therefore dependent on both flow rate and the structure and orientation of heterogeneity. In the high flow rate experiments, the equivalent relative permeability coincides with the viscous limit relative permeability, shown in Figure 1b. This is expected, since capillary pressure gradients are overcome by viscous pressure gradients at high flow rate and there are negligible capillary end effects (see the supporting information Figure S1), meaning the core-averaged saturations and voxel/slice-averaged saturations are very similar.
In Figure 1b, the core-average water saturations in the experiments range from 0.27 to 0.62, in which the relative permeability to both gas and liquid remains k rj < 0.12. This is typical of core flood experiments using low viscosity fluids where the maximum capillary pressure achievable is limited by the viscous pressure variation between the upstream and downstream ends of the rock core . The viscous limit relative permeabilities for the Bentheimer and Bunter are similar, which is due to the comparable pore structures indicated by the average capillary pressure curves for the cores, unimodal pore-size distribution (Figure 2), and average absolute permeabilities.

Capillary Pressure Heterogeneity
To characterize capillary pressure heterogeneity, we use observations of average capillary pressure characteristics and 3-D saturation maps from the core flood experiments, in an iterative numerical calibration process. The observational approach is based primarily on the work of Pini and Benson (2013a, 2013b), whereby an initial estimate of capillary pressure heterogeneity can be derived from the observed saturation heterogeneity in steady state two phase core floods. This information is used in an iterative calibration of a numerical model of the rock core. This builds on the approach of Krause (2012), Krause et al. (2013), and Krause and Benson (2015) by incorporating multiple fractional flows and total flow rates into the iterative calibration, assuming the viscous limit relative permeability is known. Distinct from past studies, we then make use of the numerical model in the derivation of equivalent relative permeability.
Small core samples are drilled from the larger rock cores after steady state core flooding for use in MICP testing to obtain capillary pressure curves which are representative of the average curve for the whole core. Figure 2 shows the MICP data which has been scaled by interfacial tension to represent the appropriate fluid pairs used in the experiments, along with fitted Brooks-Corey curves with parameters detailed in Table 1. We use the Brooks-Corey model for capillary pressure P c in this work, due to its applicability over a wide range of porous media systems (Pini & Benson, 2013a): where P e is the entry pressure, S wirr is the irreducible water saturation, and k is the pore-size distribution factor. We use the threshold pressure of the unimodal rock cores (see the pore-size distributions in the inset of Figure 2) measured in the MICP experiments as the representative entry pressure for the average system. k and S wirr are found by minimizing the fit between the MICP data and predicted P c ðS w Þ from the Brooks-Corey model with fixed P e . This gives an estimate for S wirr which is similar to that observed experimentally at the voxel level. The parameters obtained from this fitting process are shown in Table 1, with the corresponding curves shown in Figure 2.
Saturation maps from the low flow rate experiments are used as the principal information to obtain an initial estimate of the 3-D capillary heterogeneity (Egermann & Lenormand, 2005;Krause et al., 2013;Pini & Benson, 2013b;Pini et al., 2012;Reynolds & Krevor, 2015). 3-D voxel saturation (and porosity) maps are calculated by arithmetically upscaling the high resolution X-ray scan data, to improve the precision and to ensure that the voxel size is above the minimum REV needed for Darcy's law to be applicable on the continuum scale. See Zhang et al. (2000), Fernandes et al. (2012), Pini and Madonna (2016) for discussion of the REV for various sandstones. We upscale the raw scan resolution by a factor of 8 3 8 3 1 in the x-y-z directions, to achieve voxel sizes above the REV for permeability & porosity in sandstone cores (see Table 1).
As an initiating guess in the calibration process, we assume the capillary pressure in each transverse slice of the core is constant as in Krause et al. (2013) and Krause and Benson (2015). Under this assumption, the variation in capillary pressure from the average Brooks-Corey curve for each voxel i is then given by a unique scaling factor j i : where S i is the water saturation at voxel location i,P c ðS i Þ is the average characteristic capillary pressure, and P c;i ðS i Þ is the individual voxel capillary pressure. Here we effectively assume the voxel level capillary pressure scales in the core with the average capillary pressure. The scaling parameter j i is initially obtained directly from the experimental data following the approach of Benson (2013b, 2017) (see Figure 3a). We minimize the following objective function to find j i : where N v is the total number of voxels, N f is the total number of fractional flows,P c is the average capillary pressure curve, andS exp j is the slice-average saturation containing the voxel saturation S exp ij . Sðj iP c 5P c ðS exp j ÞÞ refers to the saturation predicted with the new curve j iP c at the slice-average capillary pressure and saturatioñ P c ðS exp j Þ. The minimization process is applied to all voxels, resulting in a 3-D map of the scaling factor j i for each core, shown in Figure 3b. This initial estimate already highlights the capillary layering that is present in both cores with variations Dj i % 0.5 resulting in DP e between layers of around 1-1.5 kPa, which would appear ''homogeneous'' in single phase tracer tests (Reynolds & Krevor, 2015). In the Bentheimer core, there is a high entry pressure layer parallel to the principal flow direction (z axis), more prominent at the downstream end of the rock core than in the upstream end. In the Bunter sandstone, there are repeated layers oblique to the z axis, some of which pinch out against the confining sleeve holding the rock core.
Once the initial estimate of capillary pressure heterogeneity is obtained, the assumption of a constant capillary pressure in each slice (which is not strictly valid when rP 6 ¼ 0) can be relaxed through numerical calibration, in order to derive an accurate estimate of capillary pressure heterogeneity (Krause, 2012;Krause et al., 2013). Through the numerical simulation of the low flow rate core flood experiments, we can update the capillary pressure heterogeneity to obtain a more representative distribution inside the core.
We represent the rock core numerically using an orthogonal, Cartesian mesh with cell dimensions equal to the upscaled voxels, resulting in a grid 14314339 for the Bentheimer, and 14314349 for the Bunter. The numerical core diameters are smaller than the physical cores (r 5 1.31 cm and 1.81 cm, respectively) to alleviate any scanning or upscaling artifacts near the outer circumference of the cores. The inlet superficial Darcy velocity (Q tot =A) is maintained at the same level as the experimental core by adjusting the volumetric flow rate Q tot accordingly.
The CMG IMEX TM fully implicit, isothermal immiscible multiphase porous media flow simulator is used to numerically simulate the core flood experiments. In this work, we use fictitious inlet and outlet slices to model the small gap (<1 mm) between the core faces and the end caps, where the boundary conditions can be fairly well approximated. At these locations, there is zero capillary pressure and the fluid saturation is controlled solely by gravity and the fractional flow. In the numerical model, we therefore use a very large permeability (>10D), constant core-average porosity, linear relative permeability and zero capillary pressure. Two injection wells are completed in the central 2 3 2 voxels of the inlet slice for each fluid phase to create a constant flux inlet boundary condition. Two constant pressure wells are used in the outlet slice to create a constant pressure outflow boundary condition. For an in-depth discussion on the choice of boundary conditions, see the supporting information. The rock-fluid properties in Table 1 are used as inputs to the simulator; average properties are assigned to the whole core, while the 3-D porosity and capillary pressure scaling factors are defined for each grid cell. Each fractional flow in the simulation runs until steady state; at least 5-10 pore volumes of fluid are injected and the steady state criteria discussed in section 2.1 above are observed.
To calibrate the 3-D capillary pressure heterogeneity in the core, we build on the works of Krause (2012) and Krause et al. (2013) by extending the iterative approach to multiple fractional flows and flow rates, operating solely on the capillary pressure. The calibration procedure uses as a guiding principle that a difference between observed and simulated saturation at a given location is due to an incorrect scaling of the capillary pressure curve. After each simulation iteration, we produce a capillary pressure-saturation pair at each fractional flow, which we update with the experimental saturation, i.e., ðP sim c;ij ; S sim ij Þ becomes ðP sim c;ij ; S exp ij Þ, shown in Figure 4a. These updated points are then used in a minimization process to find the new j i , in a similar fashion to the original experimental scaling in equation (6): where the slice-average capillary pressures have been replaced by the simulation capillary pressure at the voxel scale, P c;ij ðS sim ij Þ, and the saturations we minimize toward are the experimental voxel values, S exp ij . Using the above minimization, the experimental voxel saturation and simulated capillary pressure are used to update the scaling factor j i in an iterative fashion, shown graphically in Figure 4a. The minimization was performed using a single fractional flow in Krause et al. (2013) and is extended in this work to multiple fractional flows. Using multiple fractional flows reduces systematic error from calibration bias, discussed in more detail in section 3.
Upon iterative calibration of the capillary pressure curves, the simulation voxel saturations converge to the experimental values, with each iteration improving the result. The calibration is completed when the mismatch between the observations and simulation no longer improves, or improves marginally (i.e., less than 1%) with further iteration; see the supporting information for more details on the convergence of the scheme. Figure 4b highlights the resulting j i map for the two cores after six iterative calibrations, with the heterogeneity magnitude increased over the base characterization in Figure 3b. This is discussed in more detail in the results section.
Calibration of the capillary pressure heterogeneity completes the characterization effort, after which the rate invariant, intrinsic rock properties which describe isothermal, immiscible multiphase flow at the mm scale have been characterized. We now analyze the accuracy and robustness of the characterization by simulating the low flow core flood experiments with the numerical models.

Results and Discussion
In this section, the low flow rate core flood experiments are numerically simulated and the equivalent relative permeability calculated using the multiphase extension to Darcy's law, equation (3). These results give an indication of the predictive ability of the numerical models when compared to experimentally derived values. We use core-average saturations in determining equivalent relative permeabilities, and assume that phase pressures are equal at the core inlet and outlet of the experiments, i.e., P c 5 0, as in the simulations.
Saturation maps and equivalent relative permeabilities for the Bentheimer core simulated using different calibration schemes are presented in Figure 5. The transparent scatter points in the plots are points that lie outside the central 5/7th of the cores, representing voxels that are likely to be affected by discrepancies between how the simulation represents end effects and the physical experiment. Figure 5. Voxel saturation scatter plots between simulation and experiment at low flow rate (capillary limit) using the different calibration schemes. (a) Base case, using the initial experimental characterization without numerical calibration. (b) Iterative numerical calibration using fractions 1-6 to minimize equation (7). (c) Iterative numerical calibration using only fraction 6 to minimize equation (7). (d) Equivalent relative permeability of the Bentheimer core at low flow rate using the different calibration schemes. Iterative calibration results are shown after six iterations. Specific values for the fractional flows can be found in the supporting information.

10.1029/2017WR022282
The iterative calibration using all the available experimental data in Figure 5b shows the highest correlation with experimental data, with a statistically significant R 2 50:80. Since the scheme uses capillary pressure information across all the fractional flows that were experimentally tested, it accurately predicts the saturation across the whole range of experimental conditions. The spread of points has decreased significantly from the base case in Figure 5a, and there is also a less prominent tailing where the saturation is close to zero.
The iterative calibration approach with a single fractional flow (f6), analogous to that presented in Krause et al. (2013), is shown in Figure 5c. The calibration has resulted in a near perfect correlation between the simulated and experimental voxel saturations for f6, with almost all voxels lying within two standard deviations of the experimental error (dashed box around the identity line). The overall correlation considering all fractional flows has also improved from the base case, as the calibration using f6 still captures much of the underlying information associated with the capillary pressure heterogeneity. Although near perfect correlations can be obtained when calibrating to high gas fractional flows, the iterative approach does not perform as well when calibrated against an observation at low fractional flow. At low fractional flows, more voxels have saturations close to zero, meaning there is less information to update the capillary pressure curves and the iterative scheme converges very quickly to a result no better than the base case. Figure 6a highlights the improvement in saturation prediction using the iterative calibration scheme with all fractional flows. In the base case, iteration 0, the saturation heterogeneity is largely under-predicted, with clear gravity segregation occurring at the inlet of the core which is not present in the experiment. With the removal of the capillary equilibrium assumption in the iteratively calibrated results, Figure 6a (bottom), the capillary entry pressure is increased in the top half of the core, leading to a more homogeneous distribution of nitrogen at the inlet, as in the experiment.
In the iterative numerical calibration, Figure 6a (bottom), the outlet end effects present in the Bentheimer experiment (increasing S w saturation toward the outlet) have been slightly over predicted by the simulation. Using P c 5 0 at the outlet and allowing the capillary pressure of the voxels near the outlet to increase as per the experimental data has caused S w to increase above that in the experiment. Using a finite, nonzero value for P c at the outlet may help to reduce this error, however, the exact boundary condition is unknown, and cannot accurately be inferred from the experiment. We therefore choose to keep P c 5 0, and use the experimental data to infer the extent to which the data is affected by the boundary. These voxels can be seen in the simulated data as outliers in Figure 5b near simulated S N2 % 0:120:2. Upon further simulation, boundary affected points can be removed.
The predicted low capillary number equivalent relative permeability curves for each different calibration scheme can be seen in Figure 5d for the Bentheimer core. The iterative calibration schemes have improved the prediction over the base case, largely due to a better prediction of the core-average saturation, with predicted pressure drops remaining fairly similar. The iterative calibration using all fractional flows shows the best match with experiments across the tested saturation range, particularly at high water saturations.
The raising of the gas equivalent relative permeability (and slight lowering of the water permeability) in Figure 5d is caused by the parallel layering of the capillary pressure heterogeneity. The layering parallel to the principle flow direction allows the fluid phases to organize themselves more optimally into high and low entry pressure regions. This has the overall effect of boosting the equivalent relative permeability to the nonwetting phase as it flows primarily through regions of locally high N 2 saturation. This effect has also been reported before in previous work (Krause & Benson, 2015;Rabinovich et al., 2016;Reynolds & Krevor, 2015).
The iterative calibration scheme is now applied to the Bunter core, using the experimental saturation data from the 5 low flow rate fractional flows. Figure 7a shows the resulting simulation saturations after six calibration steps compared to the experimental data. The R 2 value of 0.56 highlights the significant correlation between experimental and simulated saturations, particularly at high fractional flows.
The lower fractional flows in Figure 7a show a systematic increase in simulated CO 2 saturation toward S CO2 50, compared to the experimental values. This trend is due to the capillary entry pressure which has been under-estimated at the low fractions, where S w % 1 in many voxels. Since the flow rate of CO 2 is only 0.02 mL min 21 for f1, many of the high entry pressure voxels remain at S w 5 1. The entry pressure is not scaled enough for these voxels, since the low fractional flow saturation provides little information toward the minimization process in equation (7), with the higher fractional flows biasing the result. This is highlighted in Figure 6b, whereby the low CO 2 saturations are not fully captured by the numerical model. However, the iterative calibration does show a clear improvement over the base case, with the saturation layering much more accurately captured. Incorporating further information for improved S w 5 1 voxel scaling is the feature of on-going research.
In Figure 7b, the simulated equivalent relative permeabilities for the Bunter are presented. The pressure drops and resulting relative permeabilities predicted by the simulation are within the experimental errors, shown in Figure 7b, although the average saturation is most accurate for the higher fractional flows. Both the gas and water predicted equivalent relative permeability drop below the viscous limit values, as in the experiments. A significant reduction in the gas or nonwetting phase relative permeability curves is common for perpendicular layered systems at low flow rate (Rabinovich et al., 2016;Virnovsky et al., 2004), where wetting phase relative permeabilities are also lowered, albeit not for the entire saturation range. The water relative permeability is also reduced, as with the parallel layering. In the case of perpendicular layering, the barriers to flow impact both fluid phases significantly.
The prediction of equivalent relative permeabilities at low flow rate highlights the effectiveness of the numerical models. At this point, it is worth discussing the experimental precision, and the effect this has on the resulting calibration efforts. The saturation precision at various levels of upscaling, along with pressure drop uncertainty is detailed in Table 2. The saturation uncertainty is calculated using the techniques described in Pini et al. (2012) and Pini and Madonna (2016), based on the standard deviation of CT numbers obtained upon subtraction of two dry-scan X-ray images. The voxel level uncertainty in saturation ranges from 3.2% to 4.5% for the Bunter and 1.9% to 2.6% for the Bentheimer, which is achieved through 83 upscaling of the raw scan resolution, 33 averaging of independent scans and a large z axis slice thickness used for each scan. Upon further upscaling to the slice or core level, the saturation precision is improved significantly. The pressure error in Table 2 is shown graphically in the form of relative permeability uncertainty in Figures 5d and 7b (2r Dp shown), whereby the errors lie largely within the data points.
Due to the high level of experimental precision in the saturation and pressure drop measurements, the uncertainty in the capillary pressure calibration is reduced to that resulting from systematic error. The main contributor of this error is likely the assumption of capillary pressure equilibrium which is the basis of the calibration efforts. We initially assume that the average saturation maps directly to capillary pressure, with changes in capillary entry pressure causing a corresponding change in local saturation.
The validity of this assumption can be assessed by comparing the Bentheimer and Bunter simulation results. The pressure and core-average saturations are predicted significantly better by the Bunter model, despite larger relative variance errors in the saturation and pressure data. The RMS relative error is 9.8% and 18.7% for the Bunter and Bentheimer, respectively, after six calibration iterations, with the difference largely in the pressure prediction. Since the low flow rate regime for the Bunter lies closer to the limit of capillary equilibrium than the Bentheimer (N c % 0:5 compared to N c % 1:5, respectively) the calibration of the capillary pressure curves is more effective, and a closer match can be obtained to the experiments. Lowering the experimental flow rate further would likely improve the calibration method, leading to enhanced characterization.
Next, we explore the possibility of calibrating the model based on observations obtained at one flow rate.

The Effectiveness of Single Flow Rate Characterization
Thus far, the characterization of multiphase flow properties in heterogeneous cores has been achieved using experimental data across multiple fractional flows and two flow rates. This is the most effective and robust way to characterize heterogeneous cores, since the viscous limit relative permeability and capillary pressure heterogeneity are most accurately derived in distinct flow regimes; viscous and capillary dominated respectively. However, practically, it is beneficial to only perform a single experiment, and characterize the core completely using a single flow rate data set. With this in mind, we investigate the effectiveness of characterizing the numerical models using data solely from the high flow rate experiments.
In the high flow rate experiments, N c is generally above the transition region for capillary-viscous dominated flow. However, there is still some level of heterogeneity in the 3-D saturation maps, which can be used to infer capillary pressure heterogeneity. We therefore use the saturation information from the high flow rate experiments to characterize the capillary pressure, and use the resulting numerical model to predict the low flow rate equivalent relative permeability.
We use fractional flows 6-10 and 5-8 from the high flow rate Bentheimer and Bunter experiments, respectively, to iteratively calibrate the capillary pressure, since these fractions provide the most information regarding the scaling of the capillary pressure curves, and cover a wide range of saturations. As the high flow rate experiments were also used to derive the viscous limit relative permeability, the equivalent relative permeability predicted at the high flow rate including capillary pressure heterogeneity falls on top of the viscous limit values, with whole-core voxel saturation R 2 values of 0.75 and 0.59 for the Bentheimer and Bunter, respectively, when compared to experimental saturations.
The calibrated scaling factors j i for the cores using the high flow rate and low flow rate experimental data are shown as a voxel by voxel comparison in Figure 8a. Viscous and capillary limit are used to refer to the calibration using high and low rate data respectively, to indicate the flow regimes to which the data sets lie closest. Both iterative calibration scheme results are shown after six iterations with their respective experimental data. In Figure 8a, the Bunter shows a correlated, but systematically lower j i prediction using the high flow rate experimental data. The Bentheimer shows less correlation, however, the majority of points also lie under the identify line, being under-predicted. The transparent points in the plot are the inlet and outlet 2/7th of the core, and show significantly worse correlation between the two calibration methods due to previously discussed boundary effects.
The systematic under-prediction in scaling factor is also highlighted by the density function plots in Figure  8b, where the density functions are clearly shifted to the left for the high flow rate calibrations. Despite the shift, the overall bimodal shape is maintained, with the peaks translating to the high and low capillary pressure regions in both cores. The capillary pressure scaling has therefore been retained in a relative sense, with the difference in entry pressure between the high and low regions the same, but with the overall magnitude altered.
The shift in scaling factors is likely due to the initial guess which assumes capillary pressure equilibrium, see Figure 3a. In the high flow rate experiments, N c is greater than in the low flow rate experiments, with the viscous pressure drop having more of a contributing factor to fluid distribution within the core. The error incurred in this initiating assumption appears to not be fully resolved through the iterative calibration against high capillary number observations.
To further assess the accuracy of the single, high flow rate characterization, the low flow rate experiments are simulated, with results displayed in Figure 9. There is significant scatter around the identity lines in Figures 9a and 9b, with voxel levels saturations showing significant deviation from the experimental values. Despite the voxel level saturation scatter, the core-averaged saturation and equivalent relative permeability are well approximated by the model, shown in Figures 9c and 9d. The relative success of the high flow rate Figure 8. (a) Comparison of the scaling factor j i calculated using high flow rate (viscous limit) experimental data, and low flow rate (capillary limit) data. Dark points are the central 5/7th of each core. Dashed blue and green lines show the correlation trends for the Bentheimer and Bunter, respectively. (b) Density functions of the scaling factor j i through the entire core, calibrated using high flow rate (viscous limit) and low flow rate (capillary limit) experimental data. Nonparametric normal kernel density functions are used to estimate the density function of j i . calibration in predicting core-averaged results is due to the fact that it still captures the relative difference in capillary pressure between the layers, even if the overall magnitude is not fully maintained. With this relative difference, we still predict the lowering of the relative permeability in the perpendicular layering case, and the raising of the gas relative permeability in the parallel layering case. If higher flow rates were chosen, the calibration would become increasingly worse as less information can be derived about capillary heterogeneity.
The total RMS error between simulated and experimental relative permeabilities for the low flow rate experiment is 19.4% and 29.2% for the Bunter and Bentheimer, respectively, significantly higher than the capillary limit characterization result of 9.8% and 18.7%, respectively. This highlights the increase in uncertainty incurred as observations underlying the calibration are obtained exclusively at one flow regime.
Care therefore has to be taken when choosing appropriate flow rates for the laboratory experiments when characterizing heterogeneity. An ideal situation is to perform two experiments, one at the highest achievable flow rate that still retains a low pore scale capillary number, e.g., n c < 10 26 , governed by allowable pressure drops and fluid recirculation rates, and the other at the lowest achievable flow rate, governed by measurable pressure drops, pump minimum limits, and time. This limits the number of experiments that have to be performed, and permits the greatest amount of characterization from the available data. Systematic errors are also kept to a minimum by allowing the viscous limit relative permeability to be accurately derived in the high flow rate nearing rS w 50, along with the capillary pressure calibration nearing rP c 50.
From these results, it is clear that the capillary pressure heterogeneity can be accurately, and uniquely calibrated using appropriate experimental data, leading to a rock core characterized with intrinsic, rate independent multiphase flow properties. We have extended the work of Krause et al. (2013) and Krause (2012) by incorporating multiple fractional flows and two flow rates into the characterization, and by using a macroscopic capillary number to guide the conditions at which the multiphase flow properties should be derived. Significant improvements in accuracy are gained by deriving properties in this way, and the predictive capabilities of the resulting numerical model are increased. Attention now turns to using the characterized numerical model in a predictive manner to capture equivalent relative permeabilities outside of the conditions tested in the laboratory.

Beyond Conventional Core Analysis
Long-standing hurdles in the analysis of heterogeneous rocks come from the impracticalities of purely laboratory based characterization. The equivalent properties are dependent on flow conditions, fluid properties, and the orientation of heterogeneity with respect to flow. With a numerical model of the rock core, we can now control these conditions to derive the appropriate properties as the basis for further upscaling.
We use the numerical models with relative permeability derived from the observations above the viscous limit and 3-D capillary heterogeneity calibrated with observations below the viscous limit. Core flood experiments are simulated across a range of capillary numbers, varied here by altering the flow rate. At each flow rate, we determine the equivalent relative permeability at the core scale by applying Darcy's law using the correct phase pressures in equation (3).
Numerical experiments are presented for three numerical models: (1) The original full length numerical model, with P c 5 0 in the end slices. These simulations are representative of a laboratory experiment inclusive of the impact of capillary end effects.
(2) The original numerical model with the inlet and outlet 1/7th removed and the end slices set toP c ðS w Þ. This removes the impacts of the boundary conditions, and provides the equivalent property relevant for describing flow in the reservoir system. (3) The same as (2) but with a rotated P c heterogeneity. As the experimental cores are drilled vertically or horizontally from a reservoir well, there is little control of the orientation of the heterogeneity with respect to the flow direction in the laboratory. By rotating the heterogeneity, we can estimate the relative permeability in a direction perpendicular to the experimental orientation, e.g., an estimate of the horizontal relative permeability can be made using a model calibrated by observations made from a rock cored in the vertical orientation.
The capillary pressure heterogeneity is rotated 908 anticlockwise in cubes of 14 3 14 3 14 voxels, shown in Figure 10. In this example, we only rotate the capillary pressure heterogeneity and porosity, assuming the core-average permeability and viscous limit relative permeability remain unchanged. Although the permeability would change in the physically rotated system, this example serves to illustrate the impact of capillary pressure heterogeneity orientation on the derivation of equivalent relative permeability, in which the general trends are independent of core-average permeability and viscous limit relative permeability.
The numerical experiments are conducted across the whole range of the fractional flow curve. We choose specific fractional flows using a distribution between f min f 12f min , where f min is determined by the minimum flow rate achievable in the numerical simulations. In CMG IMEX the minimum flow rate is %0.0007 mL min 21 , below the minimum possible in laboratory using 1000D Teledyne ISCO pumps (%0.001 mL min 21 ). A geometric distribution of fractional flow points allows population of the fractional curve with many points near the extremes, where a small change in fractional flow causes large changes in core saturation, allowing the full range of experimental saturations to be accessed.
Results at varying N c using the three numerical models for the Bentheimer and Bunter are shown in Figure  11. Quoted N c are based on f g 50:5 fractional flow, since this typically has the largest pressure drop in the experiments due to the competition of both fluids through the pore space and the relatively high effective viscosity of the system.
In Figure 11, the top row of plots using the original, calibrated numerical model clearly highlights the effect of outlet P c discontinuity (the capillary end effects) when compared to the simulations with the end effect eliminated shown in the second row of plots. The original cores show a large variation of equivalent relative permeability when N c < 10, with both water relative permeabilities significantly reduced relative to the viscous limit values. This is typical of a steady state relative permeability experiment with significant boundary effects, whereby the discontinuity at the outlet boundary has increased the pressure drop in the water phase, reducing the relative permeability (Krause & Benson, 2015). This effect is significantly reduced for the shortened core with continuous P c in Figures 11c and 11d, highlighting the importance of removing the boundary effects in the derivation of the equivalent relative permeability.
In the second row of plots in Figure 11, the equivalent relative permeability curves are shown with the viscous and capillary limit relative permeabilities. The capillary limit relative permeabilities were calculated using the steady state numerical upscaling technique described in Pickup and Stephen (2000). Using the calibrated numerical model in this way, the bounding relative permeability curves can be efficiently derived.
In Figures 11c and 11d, the equivalent relative permeability curves converge to the capillary limit at low flow rates and the viscous limit at high flow rates, covering four orders of magnitude in N c . The deviation from the bounding curves across a large range of flow regimes highlights the importance of equivalent relative permeabilities for modeling subsurface flow, where the flow regime changes significantly between the point of injection and the far field (Reynolds & Krevor, 2015).
The change in equivalent relative permeability from the viscous or capillary limit manifests itself differently for the Bentheimer and Bunter, depending on the orientation of the capillary pressure heterogeneity. The increase in gas permeability from the viscous limit in the parallel layered Bentheimer in both Figures 11a and 11c signifies it is a physical effect from the capillary pressure layering, as reported by previous authors (Krause & Benson, 2015;Rabinovich et al., 2016;Reynolds & Krevor, 2015). In the Bunter, we see a different effect, whereby both the gas and water relative permeability curves are lowered below the viscous limit as in Virnovsky et al. (2004) and Rabinovich et al. (2016). With combined boundary effects in Figure 11b, the curves have been more significantly reduced, even at high flow rates. This indicates that the impact of capillary heterogeneity in this domain prevails at high capillary number.
In the Bunter sandstone model, by eliminating the boundary discontinuities we significantly change the equivalent relative permeability. This suggests that obtaining equivalent relative permeabilities at low capillary number from laboratory observations alone would be significantly complicated by the presence of the boundary effects of the experiment. With the numerical models, the boundary effects may be eliminated, resulting in the derivation of an equivalent relative permeability representative of flow in the reservoir.
A further experimental constraint that can be removed is the structure of the heterogeneity, which can be reorientated in the numerical models as seen in Figure 10. The resulting relative permeabilities in Figures  11e and 11f show that the characteristic variation with capillary number has now switched for the two rock samples: the reoriented Bentheimer rock core now has equivalent relative permeabilities that are lower for both phases with decreased capillary number. This is because the layering is now perpendicular to the direction of flow. For the Bunter, the gas phase relative permeability now increases with decreasing capillary number. The reoriented layering is now perpendicular to the principal direction of flow. The features of the equivalent relative permeability are therefore driven primarily by the orientation of the heterogeneity, as opposed to the rock core itself.
At high fluid saturation in Figures 11c-11f, the relative permeability for that particular fluid converge to the same curve regardless of the heterogeneity structure or the capillary number, also seen in Virnovsky et al. (2004), e.g., Figure 6. The relative permeability of the other fluid changes significantly with capillary number at these points, varying by as much as 1.5 orders of magnitude in the results presented here. The functional dependence of the relative permeability on saturation and capillary number can therefore be anchored at the intrinsic values at high gas and water saturations, providing a clear route for potential modeling efforts, and incorporating the relative permeability functional form into simulation packages. We suggest that it is this form of drainage relative permeability that should be used in field-scale modeling.
The capillary number at which the saturation distribution and therefore equivalent relative permeability becomes homogeneous is generally reported in the range 10 N c 100 (Reynolds & Krevor, 2015;Virnovsky et al., 2004). In Figure 11, the transition to the viscous limit relative permeability does appear to occur when N c > 100 for all systems, even when boundary effects are considered, and represents an upper bound for homogeneity. By removing the boundary effects, we see the transition occur at much lower capillary numbers, around N c 5 10240 for the cores in Figures 11c and 11d, highlighting the care that has to be taken when interpreting the core flood data. The quoted capillary number is that at f g 50:5, and varies considerably across the range of fractional flows in each experiment, so much so that experiments can transition from the capillary dominated regime to a viscous dominated regime as f g is varied, see also Reynolds et al. (2018).
To assess the heterogeneity of the system, the scaling factor g5rðP e Þ=lðP e Þ can be used, proposed by Li and Benson (2015). For the Bentheimer and Bunter cores in this work, g50:16 and 0.17, respectively, when considering the central sections of each core. This level of heterogeneity is on the lower end of those considered in the work of Li and Benson (2015) where 0:2 < g < 0:5, and similar to the Berea core considered by Pini and Benson (2017), where g % 0:1. Extending this, we use a scaling factor based on an evaluation of saturation heterogeneity to describe the observations. In this way, the impact of the prevailing flow regime is incorporated. The saturation heterogeneity is represented by the scaling factor g Sw : where lðS w Þ is the core-average saturation and rðS w Þ is the standard deviation in saturation through the core. g Sw is presented for the f g 50:5 fractional flow at different capillary numbers for the numerical experiments in Figure 12. Here it is clear that after N c > 100, the gradient dg Sw =dN c reduces significantly as g Sw approaches zero and the 3-D saturation in the core becomes truly homogeneous. Only at this point, when g Sw 50 is the core saturation homogeneous, and the viscous limit relative permeability is derived, in practice at N c > 40 when there are no boundary effects. Using N c , g Sw , and g in this way, we can judge the level of heterogeneity in the core, and the likely impact on multiphase flow functions.
In the current work, we note that the variation in capillary entry pressures for the bimodel systems is only of the order 300-1,500 Pa, yet is enough to cause significant capillary number dependency in the relative permeability of the rock samples. Similarly, the heterogeneity factor g % 0:16 is low compared to simulation work by other authors (Fourar & Radilla, 2009;Li & Benson, 2015). We have found significant capillary number dependency in these rocks that would conventionally appear as homogeneous, with capillary number dependence likely to be prominent in almost all but the most homogeneous, isotropic cores, which are not representative of the subsurface. This highlights the shortcomings of conventional core analysis, whereby only the most homogeneous rocks are used to derive the viscous limit relative permeability. A relative permeability derived from such rock samples will have little predictive ability for describing flow in the reservoir system; the concept of the homogeneous rock core therefore has little relevance for a physics based approach to modeling multiphase flow in the subsurface.
The focus of the experimental workflow shifts under this new paradigm; instead of trying to most accurately find a condition invariant relative permeability, we now use the observations under optimized conditions to obtain the most accurate data-viscous limit relative permeability, capillary pressure heterogeneity-for the creation of a numerical rock core. From the numerical model, we derive equivalent relative permeabilities that incorporate rock heterogeneity in the relevant orientation and flow rate dependency, to model the effects of mm scale capillary heterogeneity at the cm scale and above. Using the characterized numerical rock core, we can undertake this workflow, producing equivalent relative permeabilities that have a functional dependence on capillary number, saturation and can be direction dependent, as demonstrated in Figure 11. Anisotropic relative permeabilities can occur as a consequence of capillary limit upscaling (Corey & Rathjens, 1956) and can be incorporated into modeling efforts in tensor form akin to the absolute permeability (Keilegavlen et al., 2012).
As well as deriving equivalent relative permeabilities, the measured j i distribution can be used to evaluate mm scale absolute permeabilities. These distributions could in turn inform small-scale geological upscaling procedures as in Pickup and Stephen (2000) and Trevisan et al. (2017a). The permeability of each voxel K i can be readily determined through the Leverett J-function relationship at JðS w 51Þ, i.e., K i 5 K / i =j 2 i /. The resulting porosity-permeability relationship at the voxel level for the Bentheimer and Bunter can be seen in Figure 13. This figure highlights the large spread of data around a typical power law fit, indicating that mm scale permeabilities cannot be readily evaluated using typical Kozeny-Carmen types of relationships (Krause, 2012).

Conclusions
In this work, we have analyzed the characterization of drainage multiphase flow properties on heterogeneous rock cores, using a combination of experiments and numerical simulation at the mm-m scale. The characterization was demonstrated using N 2 -DI water and CO 2 -brine systems at reservoir conditions on two sandstones with different principle orientations of heterogeneity, highlighting the applicability of the scheme across a wide range of experimental conditions. The fundamental experimental observations which informed the characterization were obtained from steady state coinjection core floods; the pressure drop through the core at multiple fractional flow rates and the associated 3-D saturation maps inside the core obtained at mm scale through medical X-ray imaging. Two experiments were performed for each core; an experiment at high flow rate, which allowed characterization of the viscous limit relative permeability and absolute permeability, and an experiment at low flow rate below the viscous limit, which allowed evaluation of the capillary pressure heterogeneity. Through scaling of the average capillary pressure curve in the rock core with a simple scaling factor j i , the capillary pressure heterogeneity at the mm scale was evaluated in 3-D using the experimentally observed saturation at low flow rate. Upon calibration of the numerical model against observed millimeter scale saturation heterogeneity, we accurately predict average pressure drop at the core scale and equivalent relative permeabilities. We also showed how a numerical model calibrated using a single, high flow rate data set could predict low flow rate experimental equivalent permeabilities with minor loss in accuracy, as long as some level of saturation heterogeneity prevailed in the high flow rate regime.
The calibrated numerical models were used to estimate equivalent relative permeabilities across capillary numbers ranging from 0:07 N c 400, roughly four orders of magnitude. We showed that the main features of equivalent relative permeabilities for parallel and perpendicular layered heterogeneous systems in the upscaled capillary limit were predicted as N c ! 0. As the fluid saturation approached unity, we also found that the relative permeabilities converged to the same curve, regardless of the capillary number. This observation provides a possible route for future modeling efforts, where a capillary number, saturation, and direction dependent equivalent relative permeability function could be anchored to the same curve as S j ! 1.
Removal of boundary effects inherent with experimental core floods, and restructuring of the heterogeneity in the final results section highlighted the applicability of the approach for advanced core characterization. The derived numerical model could be used without the limitations imposed by the experiment to derive meaningful equivalent relative permeabilities representative of the flow in the subsurface. These flow functions can then be used directly in modeling efforts to more accurately simulate plume evolution and capillary trapping associated with low potential flows in the subsurface.
The characterization presented in this paper represents a first step in obtaining a completely parameterized numerical model for use in reservoir characterization. Efforts are currently underway to incorporate trapping characteristics and associated capillary pressure and relative permeability hysteresis into this approach, allowing imbibition multiphase flow to be characterized in a similar manner.