Simulating Core Floods in Heterogeneous Sandstone and Carbonate Rocks

The characterization of multiphase flow properties is essential for predicting large‐scale fluid behavior in the subsurface. Insufficient representation of small‐scale heterogeneities has been identified as a major gap in conventional reservoir simulation workflows. We systematically evaluated the workflow developed by Jackson et al. (2018), https://doi.org/10.1029/2017wr022282 for use on rocks with complex porosity and capillary heterogeneities. The workflow characterizes capillary heterogeneity at the millimeter scale. The method is a numerical history match of a coreflood experiment with the 3D saturation distribution as a matching target and the capillary pressure characteristics as a fitting parameter. Coreflood experimental datasets of five rock cores with distinct heterogeneities were analyzed: two sandstones and three carbonates. The sandstones exhibit laminar heterogeneities. The carbonates have isotropic heterogeneities at a range of length scales. We found that the success of the workflow is primarily governed by the extent to which heterogeneous structures are resolved in the X‐ray imagery. The performance of the characterization workflow systematically improved with increasing characteristic length scales of heterogeneities. Using the validated models, we investigated the flow rate dependency of the upscaled relative permeability. The findings showed that the isotropic heterogeneity in the carbonate sample resulted in non‐monotonic behavior; initially the relative permeability increased, and then subsequently decreased with increasing flow rate. The work underscores the importance of capturing small‐scale heterogeneities in characterizing subsurface fluid flows, as well as the challenges in doing so.

• Integrating numerical modeling and experimental observations is essential to characterizing smallscale heterogeneity for subsurface flow • Capillary heterogeneity characterization in carbonates was only successful when key features were resolved in X-ray imagery • Isotropic capillary heterogeneity in carbonates results in a nonmonotonic rate-dependant relative permeability

Supporting Information:
Supporting Information may be found in the online version of this article.

TECHNICAL REPORTS: DATA
An approach to characterize capillary heterogeneity within rock cores has been developed which combines experimental and numerical methods (Krause et al., 2011(Krause et al., , 2013. The workflow uses a numerical history match of a coreflood experiment with the 3D saturation distribution as a matching target and the 3D capillary pressure characteristics as well as the absolute permeability distribution as a fitting parameter. In recent years, the approach has been further applied to a wide range of samples and developed to improve the observational basis and strengthen the iterative matching procedure (Berg et al., 2013;Hosseinzadeh Hejazi et al., 2019;Jackson et al., 2018;Krause & Benson, 2015;Krause et al., 2013;Kong et al., 2015;Ni et al., 2019;Pini & Benson, 2013a, 2013bReynolds et al., 2018). However, uncertainty still remains in characterizing more complex rocks typical of subsurface reservoirs across multiple fractional flows and flow rates, for example, across the full range of flow conditions that would typically be encountered.
In this study, we evaluate the approach described by Jackson et al. (2018) to characterize heterogeneity in rock cores across five rock samples selected for their range of heterogeneity types and length scales comprising two sandstones and three carbonates. For the sandstones we use the experimental coreflood data of Reynolds and Krevor (2015); Reynolds et al. (2018) and for the carbonate rocks, the data set of Manoorkar et al. (2021). The sandstones exhibit distinctly orientated planar bedding, whereas the carbonates are characterized by isotropic cementation of varying length scale. The range of rock types allows us to systematically evaluate the inversion approach against a systematic variation in rock structures. Using the validated models, we then investigate the impact of heterogeneity on flow behavior. We compare the rate dependency of relative permeability observed in sandstones and carbonates, and draw conclusions on the varying flow behavior.

Coreflood Observations
Previously acquired experimental datasets using five rock cores were studied, with experimental methods and data reported in Reynolds and Krevor (2015); Reynolds et al. (2018); Manoorkar et al. (2021). These data sets comprise observations from corefloods performed on two sandstones and three carbonate rock samples. Steady-state experiments were performed with the co-injection of nitrogen and DI water or 2 CO E and brine at high (HR) and low (LR) flow rates to obtain flow parameters in the viscous-limit (VL) and capillary-limit (CL) flow regimes, with 3D X-ray images taken throughout using a medical X-ray CT scanner.
The five rock samples cover a range of depositional facies and exhibit distinct heterogeneity types and lengthscales, Figure 1. The Bentheimer sandstone exhibits a simple porosity heterogeneity orientated as a single layer parallel to the flow direction. The Bunter sandstone exhibits characteristics of early diagenetic processes (Brook et al., 2003). As such, our sample features noticeable heterogeneity in porosity, which can be grouped into layers perpendicular to the axis of flow. The three carbonate samples, Indiana limestone, Estaillades limestone, and Edwards Brown dolomite, exhibit varying degrees of cementation resulting in distinct lengthscales of isotropic heterogeneity resolved by the CT scanner. The Indiana limestone has the smallest scale porosity heterogeneity ( E 1 mm). The Estaillades has an order of magnitude larger scale porosity heterogeneity ( E 1 cm). The Edwards Brown dolomite has a multi centimeter length low-porosity region toward the outlet of the core. Please see Supporting Information S1 for a detailed description of the geological setting of each sample and a summary of the experimental and modeling parameters (Supporting Information S1).

Numerical Modeling
The flow simulations were performed using a fully implicit, isothermal black oil fluid simulator (CMG ™IMEX), which uses the finite difference method to solve the governing equations. The simulations of each core used the rock and fluid properties given in Table S6 together with 3D porosity data obtained from medical CT images, processed following the standard method detailed in Withjack (1988). A detailed discussion of grid dimensions and representative elementary volume (REV) is provided in Section S3.
The viscous limit relative permeability curves were obtained with the high flow rate experiments. All viscous-limit relative permeabilities were assumed to be uniform throughout the rock domain and were modelled using the Chierici functional form Chierici (1984): for each voxel and saw no correlation, which confirmed this. Hence, the permeability was assumed to be constant throughout the samples and was set to the absolute permeability measured in a single phase flow experiment. For a detailed discussion, see Section S4.

History Match of the Numerical Models
The approach of Jackson et al. (2018) was used to history match the numerical models for all of the experimental datasets. With two of the carbonate rocks, the Indiana and Estaillades limestone, the approach was either unsuccessful or only partially successful. In these cases an extended approach was used. We summarize the approach of Jackson et al. (2018) and the extended approach below.

Approach of Jackson et al. (2018)
The 3D saturation maps of the rock cores obtained during the corefloods at all fractional flows were used to infer the heterogeneity in capillary pressure characteristics. An initial guess was made based on an inversion of the saturation maps, followed by an iterative process whereby capillary pressure characteristics were updated based on the comparison between simulated 3D saturation maps and the observations. The workflow is briefly described below.
For the initial guess, capillary pressure in each slice was assumed to be constant, that is, at equilibrium. The average saturation in that slice was assumed to map to the Brooks-Corey fit of the capillary pressure characteristic curve measured during routine core analysis (Equation 2). Voxel scale variation in the saturation within the slice was assumed to be caused by the capillary heterogeneity within the slice. From this a scaling factor  E was assigned to each voxel, adjusting the local capillary pressure characteristic curve, to minimize the mismatch between the slice-average and voxel-specific values: N is the total number of fractional flows and  , ( ) i c avg E S P represents the saturation of the average capillary pressure curve after it has been scaled (using the slice-average capillary pressure). Through this, a 3D map of the initial scaling factor was built. The scaling was then used to populate the capillary pressure characteristics of numerical simulations of corefloods using the CMG ™IMEX fluid simulator. The capillary heterogeneity was introduced alongside the porosity heterogeneity, the core-average characteristic capillary pressure behavior and the viscous-limit relative permeability.
The iterative calibration of the capillary pressure characteristics followed from this. After the first simulation, 3D saturation and capillary pressure maps were extracted and compared to the experimental values. A deviation from the experimental observations (both c E P and w E S ) was assumed to stem from an incorrect scaling parameter assigned to the voxel.  E was then updated, minimizing the mismatch between the experiment and simulation. The objective function becomes the following, where the slice-average c E P values have been replaced by the simulation values: represents the saturation of the average capillary pressure curve after it has been scaled using the individual simulated voxel capillary pressure rather than the experimental slice-average capillary pressure. More detail can be found in Section S4.

Extended Approaches Using Relative Permeability as a Fitting Parameter
The viscous-limit relative permeability curves used as input to the simulation are a major control on the results and ideally are measured with very little uncertainty in the observation. However, coreflood experiments are challenging and the measurements can be affected by sources of error (Berg et al., 2021;Krause & Benson, 2015). In this work the observations from the carbonate rocks were difficult to match following the conventional workflow. Hence we explored improvement in the history-matched simulation by additionally using the core-averaged input viscous-limit relative permeability as a fitting parameter, with the observed relative permeability at the low flow rate as an additional matching target. Please see Section S7 for a detailed description of this extended approach.

Results and Discussion
To evaluate the matching procedure, the voxel-scale fluid saturations were compared between the experiments and the simulations. In the initial approach, the high flow rate, viscous-limit relative permeability was input into the simulations. The evaluation of the predictive capability was based on the predicted observations at the low flow rate, which reflect the impact of the heterogeneities inverted from the experimental saturation data.

Sandstone Rocks
Following the workflow from Section 2.3.1, digital models for the two sandstones were generated. The relative permeabilities from the simulations align well with the experimental data for both samples, although the Bentheimer generally exhibits a closer match. These results demonstrate the replicability of the digital models: using identical modeling parameters in the unaltered iterative calibration workflow allowed us to regenerate the findings presented in Jackson et al. (2018). For a more comprehensive analysis, please see Section S5.

Carbonate Rocks
Numerical models incorporating capillary heterogeneity were created for the carbonates, Figure 2. The matching procedure and predictive capability using the conventional approach, where the experimental viscous-limit relative permeability is constrained, are assessed using Figure 3. Further error analysis is included in Section S6. The results using the extended approach, when the viscous-limit relative permeability is used as a fitting parameter, are presented in Section S7.

Indiana Limestone
As shown in Figure 2, the Indiana limestone displays the finest-scaled and largest range of variations in  E ( 

Estaillades Limestone
The length scale of heterogeneity for the Estaillades ( E 1 cm) is between that of the Indiana and the Edwards Brown. It is associated with the least variation in entry pressures with a range    0.5 0.2 E , Figure 2b. It exhibits a gradual distribution of  E clustered into relatively large regions within the core.
The voxel saturation plot displays a linear correlation between the experiment and the simulation ( 2 E R = 0.75) ( Figure 3c). The iterative calibration partially matches the saturation distribution captured in the CT scans, although they are systematically overestimated. This is also supported by the 3D saturation maps (Section S6). The digital model fails to predict the experimental low rate relative permeability (Figure 3d): rg E k is significantly underpredicted. The prediction of water relative permeability fits the experimental data better, but the curve is shifted to lower water saturations. This is consistent with the voxel saturation correlations. The experimental data displays a raised rg E k relative to the viscous-limit rg E k , suggesting that the heterogeneity enhances gas flow. This is not reproduced in the simulation. The offset in rg E k is not caused by some areas of Figure 3. An evaluation of the history match and predictive ability of models of carbonate rocks generated with the high flow rate relative permeability curve as input after five iterations. Left: voxel saturation correlation plot comparing the experiment and the simulation. The colors correspond to the individual fractional flows. Right: Computed relative permeabilities obtained from numerical simulations using the iteratively calibrated digital models as input. Top: Indiana, Middle: Estaillades, Bottom: Edwards Brown. VL stands for the viscous-limit relative permeabilities used as input to the simulations. Exp LR and HR refer to the low and high rate relative permeabilities respectively, obtained from the coreflood experiment. Sim LR refers to the relative permeabilities obtained from the simulation of the low rate experiment.
the core being assigned very high e E P, the range of  E values is small. The sensitivity to boundary conditions was also tested, see Section S8 for a detailed discussion.

Edwards Brown Dolomite
The initial matching procedure results in a large spread of voxel saturations, implying that the workflow is not able to match individual voxel saturations well (Figure 3e). However, the spatial distribution is reproduced-a high number of voxels plot closely to the 1:1 trendline. This is also observed in the 3D saturation maps (Section S6). The numerical model accurately predicts the capillary-limit relative permeabilities, significantly better than the other two carbonate samples. Overall, this suggests that the model captures the controlling heterogeneity within the Edwards Brown dolomite. The low rate experiment exhibited a raised rg E k compared with the viscous-limit curve, which is reproduced in the simulation.

Evaluation of the Workflow Applicability
Carbonates are associated with complex, often bimodal, pore systems, which affect the flow behavior (Cantrell & Hagerty, 1999;Prodanović et al., 2015). The Estaillades limestone in particular exhibited significant microporosity. To characterize this more closely, the MIP data was refitted using more suitable functional forms and the optimisation was repeated. The results displayed no visible improvement. A more extensive approach to model microporosity is to include variations in wirr E S within the sample, as presented in Kurotori and Pini (2021). This will be investigated in future work. Microporosity may also be characterized by a different set of relative permeabilities. To explore this, more detailed data collection is required, for example, inferring relative permeability from small plugs corresponding to distinct facies within the core. Furthermore, the permeability variation within the carbonates arising from the complex pore structure likely impacts the observed fluid behavior. However, to infer the permeability in these rocks, pore-based modeling for example pore-network modeling (Zahasky & Benson, 2018), or positron emission tomography (Zahasky et al., 2020) is required. Thus, our models do not explicitly characterize the 3D permeability distribution (Section S4). Lastly, the carbonate pore structure likely enables combined drainage and imbibition to occur during the coreflood experiments. In that case, basing the modeling effort solely on drainage behavior would not suffice. To explore this further, pore scale coreflood observations are required.
Non-uniqueness is a well-known issue in history matching (Berg et al., 2021). By using multiple datasets (low and high flow rate) and  w c E S P observations from all fractional flows (max. 16 for the Edwards Brown), the simulations were better constrained, further reducing the potential impact of non-uniqueness. The high volume of data guiding the workflow also differentiates it from previous studies (Hosseinzadeh Hejazi et al., 2019;Krause et al., 2011;Ni et al., 2019) and ensures the physical properties of the sample are replicated. Testing the fully empirical approach (Section S7) emphasized this. It led to the loss of predictive capability achieved with the initial iterative approach on the Edwards Brown. This underscores the importance of incorporating the physical mechanisms controlling the fluid behavior by integrating a large number of datasets. To evaluate the impact of grid convergence issues, a locally refined grid was implemented in the simulation of one of the carbonates. The resultant voxel saturation correlation displayed no improvement, hence this was eliminated as a major control on the results.
The success of the workflow correlates with the characteristic length scale of the heterogeneity in each sample. The Edwards Brown dolomite was the only carbonate sample with a successful model-its porosity heterogeneity was clearly resolved by the imagery. Thus, a factor likely responsible for the breakdown of the calibration is the extent to which the controlling features are resolved by the medical CT scanner. The Indiana and Estaillades limestone samples are known for their pore-scale heterogeneities (Ji et al., 2012;Lin et al., 2016;Muljadi et al., 2016;Nejad Ebrahimi et al., 2013;Pereira Nunes et al., 2016). Hence, the CT images fail to provide sufficient detail to infer the capillary heterogeneity in these two samples, which ultimately results in a poor characterization. Out of three samples with successful characterization-the Bentheimer, Bunter and Edwards Brown-the Bunter sandstone has the smallest heterogeneity. The perpendicular layers have dimensions of    30 30 3 E mm mm mm in the x, y, and z directions, respectively. Thus, the ratio of raw voxel to feature size is 13:1 in the x and y dimension and 1:1 in the z dimension. Similar

Water Resources Research
WENCK ET AL.
10.1029/2021WR030581 9 of 11 ratios have been observed as necessary in the accurate analysis of other X-ray imagery Lin et al., 2018).

Rate-Dependant Behavior of Relative Permeability
Using the successful models for the Bentheimer sandstone, Bunter sandstone and Edwards Brown dolomite, simulations were run at varying flow rates to investigate the rate dependency of relative permeability, Figure 4. The trend varies strongly between the sandstones and the carbonates. The perpendicular layering in the Bunter sandstone generally reduces the permeabilities of both phases, and as the total rate is increased, the relative permeabilities also increase. Comparatively, the Bentheimer sandstone, with parallel layering, allows for the phases to distribute more optimally, leading to a raised gas relative permeability relative to the viscous-limit curve.
The carbonate sample displays a different trend. Of key interest is that, distinct from the sandstones, the variation in relative permeability is non-monotonic with varying flow rate. From the lowest flow rate, the relative permeability initially increases in the wetting and non-wetting fluid phases before decreasing again toward the viscous-limit curve. This is related to the isotropic nature of the heterogeneity. This non monotonic behavior was hypothesized in Virnovsky et al. (2004). Additionally, the rate dependency is weaker than in the sandstones. This is because of the large size of the heterogeneity in the Edwards Brown. Ultimately, this controls the fluid behavior even in the high rate experiment. At varying rates, the shift in relative permeabilities is stronger for rw E k than rg E k . This suggests that the relative strength of capillary forces has a stronger impact on the water flow compared with the gas.

Conclusions
We have applied the capillary heterogeneity characterization workflow presented in Jackson et al. (2018) to five samples with varying degrees of heterogeneity. This allowed us to systematically evaluate the inversion approach against a systematic variation in rock structures, from relatively homogeneous sandstones to complex carbonates with subresolution porosity.
The workflow successfully characterized the heterogeneity in the two sandstones, where a good match in the experiment and simulation measurements was observed. In contrast to this, the successful application of the workflow to the carbonates was found to correlate with the size of the controlling heterogeneous features resolved by the scanner. The Indiana and Estaillades limestone with subresolution porosity led to unsuccessful models. Comparatively, the Edwards Brown with a large cemented region clearly resolved in the imagery exhibited a good agreement between the experiment and the simulation. This supports the hypothesis that the iterative calibration is dependent on the amount of controlling rock structure that can be resolved with the medical CT scanner.
To improve the predictions, we tested a fully empirical approach whereby the input viscous-limit relative permeabilities were used as fitting parameters. This approach failed to result in predictive capability for the impact of flow rate on saturation distribution or relative permeability in all three carbonates. Notably, it led to the loss of predictive capability achieved calculation and the tabulated data are presented in Sections S2 and S9, respectively. VL stands for the viscous-limit relative permeability.
with the initial iterative approach on the Edwards Brown. This finding emphasizes the importance of considering the physical mechanisms controlling the fluid behavior in the creation of a numerical model of the rock cores. The purely numerical approach to matching the experiments led to a complete loss of the predictive capabilities of the models.
For the Bentheimer, Bunter and Edwards Brown samples with successful models, we investigated the flow rate dependency of the upscaled relative permeability. We found that the isotropic heterogeneity in the Edwards Brown resulted in non-monotonic behavior; initially relative permeability was increased, and subsequently decreased with increasing flow rate.
The work both underscores the importance of capturing small-scale heterogeneities in characterizing subsurface fluid flows, as well as the challenges in doing so. Where imagery can sufficiently resolve heterogeneous features, the 3D history match procedure can be performed entirely within the realm of Darcy based reservoir simulators. For rocks where smaller scale features evidently place controls on the upscaled flow further techniques may be required such as the use of pore-network flow models (Zahasky et al., 2020).