CHARACTERIZATION OF HYSTERETIC MULTIPHASE FLOW FROM THE MM TO M SCALE IN HETEROGENEOUS ROCKS

. Incorporating mm-m scale capillary pressure heterogeneity into upscaled numerical models is key to the successful prediction of low flow potential plume migration and trapping at the field scale. Under such conditions, the upscaled, equivalent relative permeability incorporating capillary pressure heterogeneity is far from that derived conventionally at the viscous limit, dependent on the heterogeneity structure and flow rate, i.e. dependent on the capillary number. Recent work at the SCA 2017 symposium (SCA2017-022) demonstrated how equivalent functions can be obtained on heterogeneous rock cores from the subsurface under drainage conditions; going beyond traditional SCAL. Experimental observations using medical CT scanning can be combined with numerical modelling so that heterogeneous subsurface rock cores can be directly characterized and used to populate field scale reservoir models. In this work, we extend this characterization approach by incorporating imbibition cycles into the methodology. We use a Bunter sandstone core with several experimental CO2 – Brine core flood datasets at different flow rates (2x drainage, 1x imbibition and 2x trapping) to demonstrate the characterization of hysteretic multiphase flow functions in water-wet rocks. We show that mm-m scale experimental saturations and equivalent, low flow potential relative permeabilities can be predicted during drainage and imbibition, along with trapping characteristics. Equivalent imbibition relative permeabilities appear as a function of capillary number, as in the drainage cases. We also find that the form of capillary pressure function during imbibition has a large impact on the trapping characteristics, with local heterogeneity trapping reduced (or removed), if the capillary pressure drops to zero, or below at the residual saturation.


Introduction
During the buoyant migration of CO2 in a subsurface aquifer, the low flow potential coupled with small-scale rock heterogeneities can have a dramatic effect on the large-scale plume migration and trapping [1].Capillary pressure heterogeneities largely dominate the multiphase flow system, with small variations from REV to cm scale able to cause significant fluid re-distribution [2], affecting macroscopic flow functions such as the relative permeability [3].Under such conditions, the upscaled (cm-m), equivalent relative permeability incorporating capillary pressure heterogeneity is far from that derived at the viscous limit, dependent on the heterogeneity structure and flow rate, i.e. dependent on the capillary number [4].The upscaled, equivalent relative permeability can be used in large scale reservoir modelling to accurately incorporate the effects of small scale capillary pressure heterogeneity, without having to model the REV-cm scale heterogeneities explicitly.Capillary pressure heterogeneities also cause variation in trapping characteristics from the mm-m scale.Without effective inclusion of small scale trapping variations, the macroscopic trapping characteristics for use in field scale models can be significantly under or over predicted by traditional trapping models.Incorporating small scale capillary pressure heterogeneity into upscaled numerical models is therefore key to the successful prediction of both the plume migration and trapping at the field scale.Recent work [4][5][6] has shown how drainage equivalent relative permeabilities (at the cm-m scale) which incorporate mm scale capillary pressure heterogeneities can be effectively derived over a wide range of flow conditions, using a combination of experimental and numerical methods.High flow rate core flood experiments allow derivation of the viscous limit relative permeability, absolute permeability and porosity of the rock.Low flow rate experiments (near the capillary limit) allow accurate determination of the capillary pressure heterogeneity.These properties are then used in commercial reservoir simulators, such as CMG IMEX, to create a 3D numerical rock model.The numerical model allows derivation of equivalent relative permeabilities over a wide range of flow conditions, removal of boundary artifacts, and reorientation of the heterogeneity with respect to the flow direction.In this work, we build on the recent work of [4][5][6] by incorporating imbibition cycles and trapping into the relative permeability and capillary pressure characterization, allowing the hysteretic nature of the systems to be analysed.We focus on a Bunter sandstone core, also presented in [4,5] showing how the imbibition relative permeability and trapping can be characterized using experimental observations and numerical simulations.

Methodology
The characterization work presented here stems directly from that presented in SCA2017-022 [5] and in the recent paper [4].The Bunter sandstone sample was drilled from the Cleethorpes-1 geothermal borehole at 1312 m depth, with length 0.151m, diameter 0.0381m.It has a permeability of 2.2 D and average porosity of 0.24, with sequences of low permeability/high entry pressure (  ) barriers of length ~ 1-3cm running perpendicular to the flow direction (See Figure 1 for the   map.).Steady-state fractional flow core floods with medical X-ray CT scanning are used as the primary experimental observation to inform the characterization [7,8].Reservoir condition CO2 -Brine (1 molal NaCL) are used at a temperature of 53°C, pore pressure of 13.1 MPa and IFT of 34.7 mPa.s.We use two drainage core floods (high and low flow rate), an imbibition core flood at low flow rate, and 2 trapping core floods to construct the numerical model.A summary of the fractional flows and flow rates are given in Table 1.Drainage viscous limit relative permeabilities are calculated from the high flow rate drainage steady-state core flood, shown in Figure 1c.The irreducible water saturation  wirr = 0.821 is obtained from corrected mercury intrusion capillary pressure (MICP) data and the upscaled, voxel level (1.875mm x 1.875mm x 5mm resolution) saturation maps.A Chierici [9] functional form is used to fit a smooth relative permeability function through the experimental data points, with the 1D simulator SENDRA used to remove boundary artifacts.See [4] for detailed analysis of the parameter estimation.
To obtain the viscous limit imbibition relative permeability, we use data from the trapping experiments together with the land trapping model.The initial-residual (IR) plots in Figure 1a and b show the trapping characteristics of the core over the full saturation range, along with the experimental saturation error (2σs).The Land trapping model [10] relates the residual saturation   with the initial, turning point saturation   through a single, core averaged trapping coefficient, c: The trapping coefficient c is obtained through a leastsquares fit with the experimental data.In Figure 1a-b, the trapping can be seen to vary significantly through the core, with a large spread in the voxel level data, even when considering the inherent noise in the X-Ray derived saturations.This scatter (and that in the slice averaged data), is largely due to local capillary pressure barriers, trapping CO2 at a saturation higher than that at snap off (residual), discussed in more detail in the results section.
In the strongly water-wet system, the viscous limit imbibition relative permeability of the non-wetting phase (CO2) can be described using the drainage viscous limit relative permeability at the corresponding mobile CO2 saturation [10]: where   2  (  2 ) is the imbibition CO2 relative permeability, and   2  is the drainage CO2 relative permeability at mobile CO2 saturation   2  .The mobile saturation can be obtained from the Land trapping model: The viscous limit imbibition relative permeabilities are shown in Figure 1c, calculated using the core-averaged Land trapping coefficient of c = 1.3.We assume that for the water-wet system, the water relative permeability shows negligible hysteresis, and follows the drainage curve [10].To characterize the core-averaged drainage capillary pressure    (  ), we use threshold pressure corrected MICP data to fit a Brooks -Corey curve of the form: where   is the entry pressure,  is the pore-size distribution factor.In this work,   = 1.61 kPa and  = 1.43.The imbibition, core averaged capillary pressure curve is found in a similar way to the relative permeability, using the mobile CO2 saturation and an equivalent Brooks-Corey form, following the recent work of [6]: where    is obtained using equation ( 3), and   can be found using equations ( 4) and ( 5) at the point of flow reversal    where    =    : This form of imbibition capillary pressure assumes that    (  ) = 0 at the residually trapped saturation.The Brooks-Corey core-averaged drainage and imbibition capillary pressure curves for the Bunter sandstone are displayed in Figure 1d, along with a scaled imbibition curve which is scaled to have a finite capillary pressure at the residual   , i.e.    (  ) =   at the residual saturation.Experimental data [11] suggests that the macroscopic capillary pressure -saturation relationship may not equal zero at the residual, even in strongly waterwet systems, and has some finite value.We test the effect this may have on hysteretic flow using   = 1 2   , since   is not explicitly known for this system.As the final part of the characterization effort, we find the capillary pressure heterogeneity in the core.We use the same method as described in [4,6] where the coreaveraged capillary pressure -saturation relationship is scaled by a constant factor   for each voxel i, i.e.: , ( , ) =  i   ̅ ( , ) where   ̅ is the core-averaged capillary pressure.The initial scaling factor is found directly from the low flow rate drainage core floods.For each fractional flow, the slice averaged saturation is used to find the average capillary pressure using equation ( 4), from which the voxel capillary pressures in that slice are then known, assuming capillary equilibrium.Shifting the slice average saturations to the voxel saturations at each fractional flow, the new points can be used to fit a Brooks -Corey curve by varying   in a least squares approach; thus finding  i .This procedure results in a 3D map of  i and the capillary pressure heterogeneity.To remove the capillary equilibrium assumption (which is not strictly true as ∇ ≠ 0, but is needed as an initiating guess) the  i map can be updated using a simulated core flood.The 3D rock core is represented numerically in CMG IMEX using cells with dimensions equal to the upscaled voxels.The viscous limit relative permeabilities, core-averaged absolute permeability, 3D porosity map and 3D  i map are used as input to the simulation, along with corresponding flow velocities and injection periods equal to the experiment.Fictitious end slices are used to enforce the boundary conditions of the system; constant flux at the inlet, constant pressure at the outlet,   = 0 and linear relative permeabilities.Upon simulation of the experimental core floods, the voxel level   (  ) at each fractional flow are used to update  i , using the minimisation approach described in [4].After several iterations of this procedure,  i for each voxel converges, resulting in an accurate description of the capillary pressure heterogeneity, shown in Figure 2. The imbibition capillary pressure heterogeneity is assumed to follow the same form as the drainage, with the same scaling factor  i altering the core-averaged imbibition    at the voxel level [6].With the experimentally characterized, bounding drainage and imbibition   −   −   relationships, we use the Killough method [12] to scan between the curves when flow reversal occurs part way through a drainage or imbibition cycle.In the simulation, all voxels have the same viscous limit bounding relative permeability curves, with varying capillary pressure bounding curves scaled from the core average with  i .

Results
We now use the characterized numerical core to simulate the drainage, imbibition, and trapping core flood experiments described in Table 1.Firstly, the high and low flow rate drainage experiments are simulated, with results presented in Figure 3.The numerical core, with layered capillary pressure heterogeneity captures the change in relative permeability from the high flow rate (near the viscous limit) to the low flow rate (near the capillary limit) in Figure 3a.The reduction in both water and CO2 relative permeability at low flow rate is due to the capillary barriers, which increase the pressure drop through the system [13].In Figures 1b-d, the voxel saturation maps are directly compared with the experiments, showing that the saturation heterogeneity is well captured through the core at both flow rates, particularly at high CO2 fractional flows (f4+).The imbibition experiment is now simulated, using the same numerical core as in the drainage cases.The imbibition properties here are predicted entirely from the drainage characteristics and trapping experiments, and are not derived from separate imbibition steady-state core flood experiments.Figure 4 shows the results of the simulated imbibition experiment, using kr and Pc hysteresis with finite Pcr.In Figure 4a-b, the low flow rate imbibition relative permeabilities have been predicted, and show a similar level of error compared to that in the drainage case.Any error in the drainage prediction is thus carried forward to the imbibition cases here since the characteristics are predicted from the drainage case.It appears that the imbibition relative permeabilities are also equivalent functions, with the variation from the viscous limit caused in part by the capillary pressure heterogeneity and in part by the saturation state at the end of drainage.The water relative permeability in Figure 4b at low flow rate shows a similar profile to that found experimentally, however the simulation viscous limit water relative permeability has no hysteresis, i.e.    (  ) =    (  ).The raising of the water relative permeability from the drainage case appears to be due to the low flow rate regime that is dominated by capillary heterogeneity, as opposed to an inherent hysteretic feature of the underlying, viscous limit relative permeability function.The voxel and core averaged errors associated with the simulation results are shown in Table 2.All three hysteresis models have good saturation predictions, with statistically significant voxel R2 > 0.5.The form of capillary pressure hysteresis however has a more prominent effect on the pressure prediction, which is significantly over-predicted when    =    .At the low flow rate, the capillary pressure can be of the same order as the overall core-averaged pressure drop, meaning that unless the capillary pressure follows an imbibition path at lower overall Pc, the core averaged pressure drop can be over-estimated.The final part of this study considers the two trapping experiments.The experiments are simulated using the same numerical core as used previously, characterized based on the core averaged trapping model with c = 1.3.
The simulation results are presented in Figure 5, using the three hysteresis methods.With the inclusion of Pc hysteresis in Figures 5c-f, the level of trapping in the core has been significantly reduced, particularly in the case where the imbibition capillary pressure can reduce to zero at the residual in Figures 5c-d.After a long period of imbibition, if the capillary pressure curves reduce to Pc = 0, there is no Pc heterogeneity in the core and the trapping is governed completely by snap-off, i.e. kr drops to zero (this occurs in Figure 5d).With this form of capillary pressure hysteresis there is no emergence of local, capillary pressure heterogeneity trapping.In order to maintain a connected, mobile phase of CO2 (i.e.kr > 0), that exhibits local capillary heterogeneity trapping, a finite Pc is required.This causes build-up of non-wetting fluid behind local capillary barriers, akin to that which occurs during drainage.In experiment 1 the core flood time is shorter than in experiment 2, and the simulation does not reach a state of equilibrium (as it appears to do in the experiments), meaning local heterogeneity trapping occurs in all cases since Pc ≠ 0 at the end of the simulation.When the imbibition capillary pressure has a finite, non-zero value at the residual, local capillary pressure heterogeneity trapping emerges in the simulation, even at long time scales, apparent in Figure 5e-f.Localised accumulations of CO2 build up behind areas of high Pc, which are greater than would be predicted by the residual trapping model.Without this form of hysteresis, it appears that local heterogeneity trapping cannot occur over long time-scales, and that after large pore volumes of forced imbibition, the system would be completely described by the turning point saturations and the residual trapping model.The existence of local capillary heterogeneity trapping over large time scales still largely remains unknown.Recent large-scale tank experiments [14] have shown how local heterogeneity trapping can reduce over large time-scales, with the system becoming almost entirely described by a residual trapping model, implying Pc = 0 at the residual saturation.Contrary to this, most core-flooding experiments in heterogeneous rocks see some form of local capillary heterogeneity trapping, implying that Pc ≠ 0 at the residual saturation.

Conclusion
In this work, we have shown how hysteretic multiphase flow can be simulated in heterogeneous rocks based on models formulated from drainage characteristics and trapping experiments.Core averaged trapping models allowed capillary pressure and relative permeability hysteresis to be estimated, which when coupled with capillary pressure heterogeneity during both drainage and imbibition allowed prediction of imbibition relative permeabilities and mm scale saturations.This characterisation effort has demonstrated that equivalent, imbibition relative permeabilities emerge at low flow rate/capillary number, and they appear to be quantifiable in the same manner as drainage relative permeabilities in heterogeneous rocks.Simulations of trapping experiments revealed that the residual value of capillary pressure has a large effect on local heterogeneity trapping.If Pc reaches zero at the residual saturation, there is no local heterogeneity trapping at large time scales, and the saturations are exactly predicted from the trapping model.However, with finite Pc at the residual, local heterogeneity trapping emerges, even at capillary equilibrium, and results in accumulations of non-wetting phase.
This work was funded by NERC, Grant number: NE/N016173/1, and Equinor.We acknowledge PRORES and CMG for providing access to SENDRA and IMEX respectively.The authors thank Catriona A. Reynolds for sharing the experimental dataset.

Table 1 .Figure 1 .
Figure 1.(a-b) Experimental IR curves from trapping experiments (see Table 1).(c) Viscous limit drainage and imbibition relative permeabilities derived from high flow core floods and trapping experiments.(d) Drainage and imbibition capillary pressure curves, derived from MICP data and trapping experiments.

Figure 2 .
Figure 2. Characterized capillary pressure heterogeneity used in the numerical simulations.

Figure 5 .
Simulated IR plots of the trapping experiments using different hysteresis methods.(a -b) kr hysteresis only (   =    ).(c -d) kr and Pc hysteresis with zero Pcr.(e-f) kr and Pc hysteresis with finite Pcr

Table 2 .
Errors between simulated and experimental imbibition core floods.Voxel R 2 is the total statistical R 2 between all experimental and simulated voxel saturations.  ̅̅̅̅ and ∆ are core averaged values.