Rate dependency of capillary heterogeneity trapping for CO₂ storage

Water Resources Research

Department of Earth Science and Engineering, Imperial College London, London, UK CSIRO Energy, Clayton South, Victoria, Australia

Samuel J. [email protected]

Three-dimensional capillary heterogeneity trapping quantified in a Bentheimer sandstone at varying flow rates using medical X-Ray CT

Continuum analytical model of capillary heterogeneity trapping developed which agrees with experimental data

Trapping length and capillary heterogeneity trapping amount can be estimated using the analytical model with petrophysical properties

Abstract

In this paper, we experimentally quantify and analytically model rate dependent capillary heterogeneity trapping. Capillary heterogeneity trapping enhances non-wetting fluid trapping beyond pore-scale residual trapping through the isolation of non-wetting phase upstream of heterogeneities in the continuum capillary pressure characteristics. Whilst residual trapping is largely insensitive to the range of flow regimes prevalent in engineered reservoir settings, continuum theory anticipates that capillary heterogeneity trapping will be more sensitive to the balance of viscous and capillary forces that occur. We perform steady-state drainage and imbibition multiphase flow experiments at varying flow rate on a layered Bentheimer sample with in situ medical X-ray CT scanning to quantify saturation. Saturation discontinuities are observed upstream of capillary pressure barriers as a result of capillary pressure discontinuities, trapping the non-wetting phase at a saturation greater than pore-scale residual trapping alone. We confirm the flow rate dependence predicted by theory whereby the relationship between the initial and residual saturations approach a 1:1 dependence as flow rate is decreased. We develop a one-dimensional analytical model to quantify the proportion of capillary heterogeneity trapping in the system and the dimensionless trapping length scale, which agrees with the experimental data and allows for rapid estimates of trapping up to the field-scale.

Plain Language Summary

To be an effective carbon removal technology, geological carbon storage requires injected carbon dioxide (CO₂) to remain securely trapped in the subsurface. Many mechanisms trap CO₂ in the subsurface, depending on the underlying geology and the relative forces present. Of interest here are the inherent variations in rock type and structure resulting in changes in CO₂ migration and accumulation. These subsurface variations lead to increased CO₂ trapping (capillary heterogeneity trapping) compared to a uniform site, generally enhancing storage security. However, as shown in this work, the amount of capillary heterogeneity trapping depends on the rate. In this work, we experimentally quantify capillary heterogeneity trapping in a sandstone rock by performing core-flood experiments at different injection rates, mimicking the conditions in a field-project. We develop an analytical model which describes the capillary heterogeneity trapping in the system based on the rock properties, and find it agrees well with our experimental data. Our model can be used in the field to estimate the trapping length and proportion of capillary heterogeneity trapping, which would aid screening projects and operational design.

1 Introduction

Geological heterogeneity strongly impacts flow and trapping during CO₂ sequestration [Regnier2019, Kortekaas1985]. Recent advances in experimental core analysis show that both fluid migration and capillary trapping are sensitive to spatial heterogeneity on length scales (mm - cm) that are usually not resolved in field-scale flow simulations [Jackson2018, Krevor2011, Ni2019]. Understanding the impact of this small-scale heterogeneity, at industrially relevant scales, is necessary to accurately predict CO₂ saturation distributions and security of CO₂ storage projects.

Trapping of CO₂ by capillary forces is an important trapping mechanism during geological carbon storage. It occurs when water imbibes into rock partially saturated by CO₂ following injection. Within the pores of rocks ($\mu$m - mm), capillary forces snap-off and immobilise isolated ganglia of CO₂. For water wet rocks, the amount of residually trapped CO₂ on the continuum scale following imbibition depends upon its initial saturation at the beginning of imbibition [Krevor2015]. The Land model is a continuum representation of trapping frequently used to parameterise this relationship:

$$S_{nw,res}=\frac{S_{nw,init}}{1+CS_{nw,init}},$$ (1)

where the index ($nw$) refers to the non-wetting phase and $C$ is an empirical trapping constant [Jackson2020REV, Krevor2020Chp]. The trapping constant is typically evaluated from imbibition core floods (e.g. [Pentland2011], [Iglauer2011]).

Capillary heterogeneity trapping is an additional trapping mechanism resulting from heterogeneity in the continuum capillary pressure on the mm to m scale. These changes in capillary pressure result in local barriers to flow within the rock structure. They cause mm and cm scale increases in non-wetting phase saturation upstream of capillary barriers, accumulating as a continuous phase which remains effectively immobilised. This behavior can result in deviations from residual trapping models, such as the Land model [Jackson2020REV], [Harris2021]. We note that capillary heterogeneity trapping has been called a variety of different terms, including local capillary trapping and capillary pinning [Saadatpoor2010, zhang2025]. We use the term capillary heterogeneity trapping here, since it encapsulates the mechanism by which the trapping manifests; heterogeneity in capillary pressure characteristics.

Capillary heterogeneity has the potential to increase the security and capacity of CO₂ storage in underground aquifers [Krevor2011, Harris2021], however the evolution of this trapping with flow rate is poorly understood [Reynolds2018, Niu2015]. Continuum theory predicts a decrease trapping with increasing flow rate due to the relative increase in viscous force over capillary forces [Debbabi2017, Ringrose1993]. In part because of this sensitivity, capillary heterogeneity could explain deviations in measured initial-residual curves observed in the literature [Dance2016, Bachu2013, Al-Menhali2016, Krevor2012] and scatter in the trapping data that is often attributed to experimental error [Krevor2015, Reynolds2018, Ni2019].

The aim of this paper is to investigate the combined effects of continuum-scale capillary heterogeneity and flow rate on core flood evaluations of trapping. A series of drainage and imbibition core-flood experiments at three different flow rates, imaged using medical X-ray CT, demonstrate the rate-dependency of measured initial-residual relationships. An analytical model is derived to analyse these behaviours, highlighting the key parameters to include when modeling and upscaling capillary heterogeneity trapping.

2 Analytical Model of Capillary Heterogeneity Trapping

The 1D semi-analytical model presented below is developed from the conventional continuum theory of multiphase flow in porous media to describe the capillary heterogeneity trapped saturation distribution during steady state imbibition through a layered system. It is similar to the models presented by Duijn et al. (1995) and Dale and Ekrann (1997) that describe the saturation profile within 1D heterogeneous oil reservoirs [Dale1997, Duijn1995] but extends them to the trapping stage of a CO₂ storage project where capillary, viscous and gravity effects may be important. The model is used to parameterise saturation discontinuities observed upstream of capillary pressure barriers, ensuring the capillary pressure continuity condition is fulfilled [Debbabi2017]. We present an overview of the derivation here, referring the reader to the supporting information for further details.

The model considers a 1D medium formed of 2 homogeneous layers with flow perpendicular to those layers (Figure 1). Each layer has different properties, with the layer closest to the injection point having the lower capillary pressure. It is representative of an element in a composite core and vertical flow through two layers in a stratified aquifer [Dale1997].

Carbon dioxide already exists in the system at the start of imbibition, unlike in the primary drainage case. During imbibition, water displaces the pre-existing CO₂. Some of this CO₂ becomes trapped upstream of the heterogeneity interface due to the contrast in capillary pressure between the layers. Eventually, at steady state, only water is flowing. At this time the CO₂ in the upstream layer is immobilized. This is capillary heterogeneity trapping.

Figure 1: (a) Diagram of the model system with length L. A region of high capillary entry pressure (layer 2) lies above a region of low capillary entry pressure (layer 1), with heterogeneity length h. (b) The Brooks-Corey capillary pressure - saturation functions for the two layers of different entry pressure are shown. During the steady state imbibition of water, the critical CO₂ saturation (S_critical) forms in the low capillary entry pressure layer, upstream of a capillary pressure barrier, as a consequence of capillary pressure continuity.

Using fractional flow theory, the saturation gradient in a 1D homogeneous porous medium at steady state when only water is flowing, is given by:

$$\frac{dS_{w}}{dx}=\frac{\frac{q_{T}\mu_{w}}{k_{rw}K}+\Delta\rho g}{P_{c}^{\prime}}$$ (2)

where $q_{T}$ is the total wetting phase Darcy flow velocity, K is the absolute permeability, $k_{rw}$ is the wetting-phase relative permeability, $\mu_{w}$ is the wetting-phase viscosity, $\Delta\rho$ is the density difference between the fluids, g is the component of gravitational acceleration in the direction of flow, and $P_{c}^{\prime}$ is the gradient in capillary pressure with respect to wetting-phase saturation $(\frac{dP_{c}}{dS_{w}})$. The capillary pressure is defined as, P_c(S_w) = P_CO2 - P_w.

We can obtain an equation describing the distribution of CO₂ trapped beneath the interface between the layers by integrating Equation 2. This gives an expression for x, the distance of the saturation from the heterogeneity, as a function of $S_{w}$ in the layer closest to the injection point.

$$x=h-\int_{S_{w}}^{S_{w,critical}}\frac{P_{c1}^{\prime}}{\frac{q_{T}\mu_{w}}{k_{rw}K_{1}}+\Delta\rho g}ds_{w}$$ (3)

where $S_{nw,critical}$ = 1 - $S_{w,critical}$, is the critical CO₂ saturation immobilized at the interface between layer 1 and layer 2, following water imbibition. The saturation spatial distribution results in a capillary pressure gradient which ensures flux continuity. The extent of the trapped saturation behind the capillary barrier depends on the functions of the upstream region, in which the trapping takes place.

The critical CO₂ saturation trapped at the interface ($S_{nw,critical}$) is determined by applying the capillary pressure continuity boundary condition (see previous work by [Dale1997, Duijn1995, VanDuijnDeNeefM.1996, VanDuijn2016]). In order to do this we need a functional form for the imbibition capillary pressure in each layer. There is a limited literature on functional forms for imbibition capillary pressure so, for simplicity, we use the Brooks-Corey model (originally derived for drainage) to parameterize the imbibition capillary pressure relationship, $P_{c}=P_{e}S^{*}_{w}{}^{-1/\lambda}$ where P_e is the capillary entry pressure and $\lambda$ is the Brooks-Corey parameter. $\lambda$ depends upon pore-size distribution and is a measure of pore-scale heterogeneity within a system [Ni2019, BC]. The variable S_w* is the water saturation normalised with respect to the irreducible water saturation S_wirr, $S_{w}^{*}=\frac{S_{w}-S_{wirr}}{1-S_{wirr}}$. Experimental data suggest that the shape of imbibition capillary pressure curves is not significantly different from drainage curves [Kleppe1997], [Wang2016], although the analysis is valid for other functional forms.

Using this we can calculate the capillary pressure at the interface between layers 1 and 2,

$$P_{c}^{threshold}=P_{e2}(S^{*}_{w2,x=h})^{-1/\lambda}$$ (4)

$$P_{c}^{threshold}=P_{e1}(S^{*}_{w1,critical})^{-1/\lambda}$$ (5)

Hence the critical wetting phase saturation at that interface is calculated

$$S^{*}_{w1,critical}=(\frac{P_{e2}}{P_{e1}})^{-\lambda}S^{*}_{w2,x=h}$$ (6)

The maximum capillary heterogeneity trapped saturation thus depends on the ratio between the upstream and downstream capillary entry pressures. The CO₂ saturation in the upper layer (layer 2) is assumed to be displaced to its residual value on imbibition in the simple two-layer system ($S^{*}_{w2,x=h}=S^{*}_{w,res2}$).

From Equation 3 we see that the extent of the trapped saturation behind the capillary barrier also depends on the wetting phase relative permeability function ($k_{rw}$). In this work we assume universal relative permeability curves over the domain, a common modelling assumption supported by Burdine’s theory [Li2015]. They are parameterised using the Corey power law model, $k_{rw}=k_{rw}^{e}S^{*}_{w}{}^{n}$, where n is the Corey parameter and k_rw^e is the end point wetting phase relative permeability [Niu2015]. Substituting this functional form, along with the Brook’s Corey model for capillary pressure, into Equation 11 gives,

$$x=h-\int_{S_{w}}^{S_{w1,critical}}\frac{\frac{P_{e1}}{\lambda}s^{*}_{w}{}^{-1/\lambda-1}}{\frac{q_{T}\mu_{w}}{K_{1}k_{rw}^{e}s^{*}_{w}{}^{n}}+\Delta\rho g}d{s_{w}}$$ (7)

Equation 7 is solved numerically, although it can be solved analytically if simple linear functions for capillary pressure and relative permeability are chosen (that is, if n and $\lambda$ are 1).

2.1 Dimensionless trapping length-scale

Equation 7 provides the system parameters which influence capillary heterogeneity trapping and so can be used to evaluate under which conditions laboratory measurements of initial-residual trapping relationships may be affected by capillary heterogeneity within the core. Reducing these parameters to dimensionless numbers provides insights into the relative importance of gravity, capillary and viscous effects on the core flood and hence on the initial-residual trapping measurements. This indicates at which storage sites and conditions capillary heterogeneity trapping is expected to be significant.

Traditional dimensionless numbers often fail to describe systems displaying capillary heterogeneity trapping due to the exclusion of key parameters incorporating heterogeneity [Ni2021]. Using the analytical model, we define a novel dimensionless number, the dimensionless trapping length scale,

$$\overline{x_{T}}=\frac{x_{T}}{h}=\frac{h-x_{S=S_{res1}}}{h}$$ (8)

where h is the thickness of the high permeability layer in which the trapping takes place (layer 1), x_T is the distance over which capillary heterogeneity trapping occurs, and $x_{S=Sres1}$ is the location at which the capillary heterogeneity trapped gas saturation becomes zero (i.e. the saturation equals the residual saturation).

Defining the viscous to capillary number $N_{v/c}$ as

$$N_{v/c}=\frac{\lambda q_{T}\mu_{w}}{K_{1}k_{rw}^{e}P_{e1}}$$ (9)

and the macroscopic gravity to capillary number $N_{g/c}$ as

$$N_{g/c}=\frac{\lambda\Delta\rho g}{P_{e1}}$$ (10)

Substituting for these into Equation 7 gives

$$\overline{x_{T}}=\frac{1}{h}\int_{S_{res1}}^{S_{w1,critical}}\frac{s^{*}_{w}{}^{-1/\lambda-1}}{\frac{N_{v/c}}{s^{*}_{w}{}^{n}}+N_{g/c}}ds_{w}.$$ (11)

The dimensionless trapping length incorporates key parameters such as heterogeneity contrast and capillary pressure characteristics. This equation shows that the larger the density difference then the larger $N_{g/c}$ and the smaller the dimensionless trapping length. If the flow is horizontal then the trapping length depends primarily on the viscous to capillary number: the larger the flow rate then the smaller the trapping length.

3 Experimental

We investigated the impact of capillary heterogeneity on initial-residual trapping in a real system by performing a series of core floods in an heterogeneous Bentheimer core. A medical X-ray CT scanner was used to determine the porosity distribution and capillary heterogeneity trapped gas saturation along the core. Additional experimental details and figures are provided in the supporting information.

3.1 Materials

The Bentheimer sandstone core was chosen as it had layers perpendicular to the axis of the core that were clearly visible to the naked eye. The average core porosity was $0.214\pm 0.003$, determined from the CT numbers using the method described in Krevor et al. (2012) [Krevor2012]. Variations in porosity along core length are shown in Figure 2. In the greyscale image, bright regions highlight the major heterogeneities, probably due to iron oxide deposits. The average absolute permeability ($K$) was estimated to be 0.94 D [Reynoldsthesis2016].

Figure 2: Grey scale visualisation of CT number for a central YZ slice (top) and slice average porosity along the normalised core length (bottom), calculated from medical-CT scans of the Bentheimer core.

N₂ and DI water were used as working fluids in analogy to a supercritical CO₂-brine system [Niu2015]. The water viscosity, nitrogen viscosity and interfacial tension were estimated respectively as $\mu_{w}$ = $1.06\times 10^{-3}$ Pa$\cdot$s, $\mu_{N_{2}}$ = $1.96\times 10^{-5}$ Pa$\cdot$s and $\gamma_{N_{2}-w}$ = $0.067$ N/m for the experimental conditions[Jackson2020].

3.2 Method

N₂-DI water core-floods were performed to measure the initial-residual saturation relationship at elevated pressure (10 MPa pore-pressure) and ambient temperature. First a drainage core-flood was performed where 100% N₂ was injected. The fluid distribution was then imaged once the system reached steady state. Due to the capillary end effect and the impact of heterogeneity, a range of initial non-wetting phase saturations were observed within the core [Niu2015]. An imbibition core-flood with 100% water injection was then carried out to obtain the residual saturation. By correlating the final saturation in each voxel to the saturation prior to imbibition, the initial-residual saturation relationship was obtained [Krevor2015].

This experimental sequence was performed at three imbibition flow rates with a constant drainage flow rate between each run, Table 1, to isolate the impact of imbibition rate on capillary trapping. The corresponding pore-scale capillary number, $N_{c}=\frac{q\mu}{\gamma}$, where $q$ Darcy velocity and $\mu$ viscosity of the invading phase was also calculated for each experiment. Full details of the experimental methodology are provided in the supporting information.

	Rate [cm3min-1]		Capillary number [-]
Run	Drainage	Imbibition	Drainage	Imbibition
1	10	0.05	4.1$\times$10-8	1.$\times$10-8
2	10	0.1	4.1$\times$10-8	2.2$\times$10-8
3	10	0.5	4.1$\times$10-8	1.1$\times$10-7

Table 1: Injection rates used in the initial-residual trapping experiments and associated experimental pore-scale capillary number.

The saturation distributions were measured using medical CT scanning with voxel size 0.09 mm (x,y) and 1 mm (z - slice thickness). In order to reduce uncertainties in voxel-scale saturations, the voxels were coarsened by averaging over 20 voxels in the x-y directions. The resulting coarsened voxels had dimensions of 2 $\times$ 2 $\times$ 1 mm. These values were chosen to give the coarsened voxels well-defined Darcy-scale flow properties, with representative elementary volume (REV) similar to those in the literature [Jackson2020REV, Herring2013, Pini2013]. The voxel-scale N₂ saturation values reported in this study are from these coarsened voxels.

Jackson et al. (2020) also carried out initial-residual trapping core-flood experiments on a heterogeneous Bentheimer sample similar to the core used in these experiments [Jackson2020REV]. Their cores had similar perpendicular layered heterogeneity as they came from the same slab. Therefore, the results of Jackson et al. (2020) are used as a literature comparison to this study.

3.3 Quantifying Capillary Heterogeneity Trapping

The 1D analytical model outlined was applied to the results of the core-flood experiment to test its utility in real 3D systems. Equation 7 was calculated with the gravity term excluded as the core-flood was performed in the horizontal orientation. The input parameters for the model, detailed in Table 2, were taken from the experiment and from Jackson et al. (2020) [Jackson2020REV].

Parameter	Value	Unit
$K$	9.87x10-13	m2
$k_{rw}^{e}$	1	-
$n$	4	-
$\mu_{w}$	0.001	Pa$\cdot$s
$\lambda$	2.7	-
$S_{wirr}$	0.08	-
$C$	1.24	-
$h$	0.07	m
$q$	7$\times$10-7 - 7$\times$10-6	ms-1
$P_{e}$	1.2 - 3.7	kPa
$\Delta P_{e,h}$	0.5 - 1.5	kPa

Table 2: Parameters used in the analytical model, with values obtained from the experiment and from Jackson et al. (2020) [Jackson2020REV]. $P_{e}$ is the capillary entry pressure of the layers and $\Delta P_{e,h}$ is the change in capillary entry pressure over the heterogeneity.

The capillary entry pressure of the different layers in the core was not directly measured. Instead, based off Brooks-Corey predictions from Jackson et al. (2020), the capillary entry pressure is estimated as 3.7 kPa and 1.2 kPa for drainage and imbibition respectively. In addition, the change in the capillary pressure over the heterogeneity was estimated as 1.5 kPa [Jackson2020REV]. Due to the uncertainty in assigning capillary entry pressure to the different layers within the core, analytical solutions were calculated for a range of values.

To understand the factors influencing capillary heterogeneity trapping, the total trapped saturation was compared to that trapped by pore-scale residual saturation alone, with differences attributed to capillary heterogeneity trapping. The pore-scale residual trapping was estimated from Equation 1 using $C_{max}$, the maximum Land trapping coefficient obtained experimentally. The proportion of trapping resulting from capillary heterogeneities was then quantified as,

$$H\%=\frac{\textrm{System average capillary heterogeneity trapped saturation}}{\textrm{System average total trapped saturation}}\times 100.$$ (12)

4 Results

4.1 Initial-residual trapping relationship

Examination of Figure 2 shows that there is a zone of reduced porosity approximately 0.55 of the distance between inlet and outlet. Capillary heterogeneity trapping is observed experimentally within the heterogeneous Bentheimer core, with an increase in non-wetting phase saturation upstream of this location. Figure 3 plots the voxel-scale Land trapping coefficient in 3D for the three different flow rate experiments. A clear rate dependency in the trapping relationship is observed. These results show that as flow rate is decreased, the Land trapping coefficient decreases upstream of the heterogeneity as a result of capillary heterogeneity trapping. The capillary heterogeneity trapped saturation remains stable for many pore volumes of brine injected, establishing itself as a long term trapping mechanism rather than a transient effect.

Figure 3: Voxel-scale Land trapping coefficient within the Bentheimer sample for three flow rate experiments, displayed in 3D. Flow is from left to right in all images.

The experimentally measured initial-residual trapping relationships are a function of water imbibition rate, as expected from the analytical model. Figure 4a plots the initial-residual N₂ saturation for three different imbibition flow rates, with the Land trapping relationship defined according to Equation 1. The maximum Land trapping relationship, $C=1.24\pm 0.13$ with one standard deviation uncertainty, best characterises the primarily pore-scale phenomenon of residual trapping [Reynolds2018]. This is calculated from the high flow rate (0.5 cm³min^-1) experiment, which obeys a Land trapping relationship with minimal scatter, implying the results are a good descriptor of the pore-scale residual trapping relationship. Further plots of the voxel level initial-residual relationship are given in the supporting information.

Deviations in the residual saturation above the maximum Land trapping relationship are attributed to capillary heterogeneity trapping. Figure 4a shows that as the imbibition flow rate is decreased, there is a bigger difference between the remaining saturation and the maximum Land trapping relationship, due to the increase in trapping by capillary heterogeneities.

Figure 4: (a) Slice average initial-residual N₂ saturation for three flow rate experiments in the Bentheimer sample (blue, orange and yellow markers). The maximum Land trapping relationship, calculated as the slice average relationship from the high flow rate (0.5 cm³min^-1) experiment, is plotted as a solid line ($C=1.24\pm 0.13$). In addition literature comparisons from Jackson et al. (2020) for a homogeneous and heterogeneous Bentheimer [Jackson2020REV] are shown (grey and black markers). The average trapping relationships for the homogeneous $C=1.34$ and heterogeneous $C=0.78$ literature samples are also shown (dotted lines). (b) Slice average initial-residual N₂ saturation for the imbibition rate 0.05 cm³min^-1 experiment. The colour gradient shows where the saturations occur along the core, from inlet ($x$=0) to outlet ($x$=$L$). The Land trapping relationship upstream of the heterogeneity $C=0.55$ (dashed line), as well as the maximum $C=1.24$ (solid line) and average $C=0.86$ (dotted line) Land trapping relationships are displayed

The experimental results were compared to the trapping experiments carried out on a homogeneous and heterogeneous Bentheimer sample in Jackson et al. (2020) [Jackson2020REV]. Figure 4a shows the high flow rate (0.5 cm³min^-1) experiment displays a similar trapping behaviour to the Jackson et al. (2020) homogeneous sample. The Land trapping relationship from the high flow rate experiment is $C=1.24\pm 0.13$, whilst the homogeneous experiment from Jackson et al. (2020) has a best fit Land correlation (slice average) of $C=1.34\pm 0.12$, with both curves falling within the 1 standard deviation of the other. In contrast the low flow rate (0.05 cm³min^-1) experiment displays a similar trapping behaviour to the Jackson et al. (2020) heterogeneous experiment. Jackson et al. (2020) used a lower Land trapping coefficient of $C=0.5$ to parameterise the trapping in this sample. This trapping relationship is comparable to our low flow rate experiment (0.05 cm³min^-1) fitted upstream of the heterogeneity, $C=0.55\pm 0.05$.

Comparing the results of these experiments to the analogous homogeneous and heterogeneous literature samples highlights the variable impact of heterogeneity on trapping. The changing nature of capillary heterogeneity trapping in our experiment suggests that, depending on the flow conditions used, the initial-residual trapping relationship may behave more similarly to a heterogeneous or homogeneous sample. This highlights the difficulties in measuring trapping relationships in the laboratory and then applying them to field-scale projects where imbibition flow rates may be different.

At low flow rates, the contrast between the maximum and minimum Land trapping coefficients determined from different sections of a heterogeneous core is much higher. Figure 4b highlights the difference between extracting the maximum Land trapping coefficient ($C=1.24\pm 0.13$) and average Land trapping coefficient ($C=0.86\pm 0.37$) from experimental data at low flow rate (0.05 cm³min^-1). A clear variation in the trapping relationship is observed for initial saturations greater than approximately $S_{N2}=0.5$, corresponding to regions upstream of the heterogeneity. The average Land trapping coefficient for the region upstream of the heterogeneity is $C=0.55\pm 0.05$, with the higher remaining saturation indicative of capillary heterogeneity trapping. Isolating this region, the small uncertainty in the trapping relationship (9%), suggests this data region may be used to capture the capillary heterogeneity trapping relationship. Examining the data nearest the core inlet in Figure 4b suggests a linear trapping relationship may be better suited to parameterise this region [Krevor2011].

Based on these results, we suggest that a high standard deviation in the Land trapping coefficient determined from low flow rate experiments may be indicative of capillary heterogeneity trapping within the core. For example in our results, one standard deviation in the average Land trapping coefficient is 44% for low flow rate experiment (0.05 cm³min^-1), compared to 10% for the high flow rate experiment (0.5 cm³min^-1). This suggests that the degree of uncertainty in the average trapping relationship may be used as a measure of the potential presence of capillary heterogeneity trapping in a system.

4.2 Rate dependent capillary heterogeneity trapping

The capillary heterogeneity trapped saturation is quantified by comparing the remaining saturation trapped to that trapped from pore-scale residual trapping alone. Figure 5a shows the proportion of capillary heterogeneity trapping $(H)$ within the experiments as a function of water imbibition rate. A significant uncertainty in the proportion of capillary heterogeneity trapping within the system results from uncertainty in the maximum Land trapping coefficient extracted. The uncertainty in fitting the maximum initial-residual saturation relationship (10%) is much larger than the slice-average saturation uncertainty due to experimental precision ($<0.15\%$) [Pini2012, Jackson2018].

Figure 5: (a) Proportion of capillary heterogeneity trapping ($H$%) over different flow rates, estimated for the Bentheimer experiments from the analytical model for a range of capillary entry pressures. In addition, the experimental results (The proportion of capillary heterogeneity trapping in the system as a function of water imbibition rate, with $\pm$ 1 standard deviation of uncertainty in $C$ shown.) and associated uncertainties are shown. The capillary pressure parameters in the analytical model are optimised to provide a good fit to the experimental data. (b) Analytical solution of trapped saturation compared to experimental results, for flow rate 0.05 cm³min^-1. The capillary pressure parameters in the analytical model are optimised to provide a good fit to the proportion of capillary heterogeneity trapping from the experimental data.

A good agreement is observed between the proportion of capillary heterogeneity trapping estimated from the analytical model and obtained from the experiments as a function of rate, Figure 5a. The analytical model was run over the range of rates and capillary pressures outlined in Table 2, resulting in a range in the amount of capillary heterogeneity trapping. The initial saturation in the analytical model was assumed to be the core average initial saturation ($S_{init}=0.5$). The results of the analytical model fall mainly within the uncertainty bounds in estimating the proportion of capillary heterogeneity trapping from the experimental results. This shows that the 1D model is able to estimate the proportion of capillary heterogeneity trapping within a real 3D core when the heterogeneity pattern is approximately 1D. The difference between the analytical model and experimental results is highest at high flow rates due to differences in how the pore scale residually trapped gas is handled at the boundary. Figure 5a shows the best fit between the analytical model and experimental results. This was obtained by selecting a capillary entry pressure and capillary pressure contrast between the low and high capillary pressure layers that best matched the experimental results; $P_{e1}=$ 3 kPa, $\Delta P_{c}=$ 0.8 kPa.

Figure 5b overlays the slice average saturation predicted from the best fit 1D model with the experimental results at low flow rate (0.05 cm³min^-1). In this instance, the slice average initial saturation from the experiment was used. It should be noted that whilst the 1D model provides a good match to the proportion of capillary heterogeneity trapping in the system, the model can only estimate the 1D saturation distribution because in reality the 3D core is more complex, with other small-scale heterogeneity present and a gradation in grain size at the heterogeneity.

Equation 8 defined the dimensionless trapping length as a means of predicting the proportion of capillary heterogeneity trapping in a system. The dimensionless trapping length was estimated for the Bentheimer system using the best fit capillary pressure described in the previous paragraph, Table 3. Figure 6 shows the dimensionless trapping length can estimate the proportion of capillary heterogeneity trapping within the real rock core-flood at the experimental conditions used, as well as for different rates, Brooks-Corey parameter ($\lambda$) and Corey parameter ($n$). This suggests that the dimensionless trapping length is a good measure for estimating capillary heterogeneity trapping within a system.

Rate [cm3min-1]	$x_{T}/h$	$H$% predicted	$H$% experiment
0.05	1.01	15	15
0.1	0.51	7.6	7.6
0.5	0.10	1.1	0

Table 3: Dimensionless trapping length calculated for the different experimental flow rates. The proportion of capillary heterogeneity trapping ($H$%) is predicted using Figure 6 and is shown to agree well with the proportion calculated from the experimental results, within the uncertainty in characterising the pore-scale residual trapping relationship from experiments.

Figure 6: Proportion of capillary heterogeneity trapping ($H$%) as a function of dimensionless trapping length ($x_{T}/h$). The analytical model is run over many orders of magnitude, including varying the Brooks-Corey ($\lambda$ {0.5-1.5}) and Corey parameters (n {4-7}). The dimensionless trapping length is estimated for the Bentheimer experimental system with optimised capillary pressure, plotted against the proportion of capillary heterogeneity trapping calculated from the experimental results. Error bars indicate uncertainty in the proportion of capillary heterogeneity trapping calculated from experiments, resulting from uncertainty in defining the pore-scale residual trapping relationship.

5 Discussion

5.1 Implications for core-flood characterisation

The experimental results show that capillary heterogeneity trapping depends upon the flow rate used in the core flood: the higher the flow rate then the lower the amount of capillary heterogeneity trapping. This is consistent with the 1D model for layered heterogeneous systems (Equation 7). In contrast, in homogeneous rocks, Niu et al. (2015) showed that the average Land trapping coefficient is rate independent [Niu2015]. Thus, a key indicator that an experiment is influenced by capillary heterogeneity is that the initial-residual trapping curves change with flow rate. The 1D model (Equation 11) shows that, for horizontal corefloods, the ratio of capillary to viscous forces is key to determining the amount of capillary heterogeneity trapping. In addition, the critical saturation trapped at the heterogeneity boundary is controlled by the ratio of the capillary entry pressures of the layers.

The average trapping relationship is typically extracted from total initial-residual core-flood measurements [K43, Krevor2012, Dance2016]. This core-scale trapping coefficient, in most cases, is an upscaled value that incorporates capillary heterogeneity trapping within the core. Jackson et al. (2019) have shown that a Land trapping model determined in this way may over predict the disconnected pore-scale residually trapped saturation in heterogeneous cores, whilst under-predicting the connected saturation caused by high entry pressure layers building-up the non-wetting phase upstream [Jackson2020REV].

Given that most rock cores probably contain small-scale heterogeneity, care is needed when undertaking core measurements of the initial-residual trapping curve. One approach is to perform experiments at a high capillary number to suppress capillary heterogeneity trapping. Alternatively, based on the findings of this work, the maximum Land trapping coefficient (rather than the average) may be extracted from experimental data to capture the impact of pore-scale residual trapping alone. Reynolds et al. (2018) proposed a similar approach of characterising a lower bound on the trapping, associating trapped saturations greater than the Land model to local capillary heterogeneities [Reynolds2018]. Characterising capillary heterogeneity trapping is challenging due to the dependence on system conditions. Experiments should be carried out under characteristic conditions to quantify the impact of heterogeneity on trapping in the field. Alternatively, initial-residual trapping experiments carried out over different flow rates demonstrate the range in potential capillary heterogeneity trapping.

5.2 Implications at the field-scale

The purpose of measuring initial-residual trapping relationships in the laboratory is to provide input data for field-scale simulations. Such simulations are used to design storage schemes to maximise storage security, including evaluating the amount of residually and capillary heterogeneity trapped CO₂. In practice imbibition does not occur at a constant rate during field scale CO₂ storage as the imbibing water is driven both by capillary forces and buoyancy [Mouche2010]. The relative magnitude of these varies across the field and with time, depending on the initial CO₂ saturations resulting from injection and the rate of CO₂ plume segregation from the water post injection. Using laboratory measurements of average initial-residual curves determined from low rate experiments may result in over-prediction of residually trapped CO₂ saturations. Conversely, using initial-residual curves determined at high flows will mean the simulation may underestimate the amount of core-scale capillary heterogeneity trapped CO₂.

A further consideration is that capillary heterogeneities occur in aquifers from all scales from the sub-core scale to greater than the metre scale [Jackson2020]. These are not typically resolved in field scale simulation where the vertical dimensions of grid blocks may be 1-10m. These sub-grid block and greater than core scale heterogeneities will result in further capillary heterogeneity trapping. Upscaling is needed to convert laboratory measurements of initial-residual trapping data to the field. Given the rate dependency of capillary heterogeneity trapping we hypothesize that these upscaled functions may need to incorporate a rate dependent Land trapping parameter. This could be achieved by proving tables of the Land trapping coefficient for different rates or possibly devising a different formulation altogether.

Capillary heterogeneity trapping has important consequences in field-scale projects, particularly at sites without a laterally extensive caprock. In these composite confining systems, capillary entry pressure variations are used to baffle flow, controlling the overall footprint of the CO₂ plume [Bump2023]. Therefore, identifying the mechanism by which saturation is trapped at the core-scale and extracting the appropriate trapping relationship when upscaling to the field is critical. For example, considering an initial saturation in the range $S_{CO2,init}$ 0.3–0.6, extracting the average or maximum Land trapping coefficient from the experimental results would correspond to a difference of 9–15% in the capillary trapped saturation. In a site injecting 1 Mt/yr this would equate to 90–150 kt of CO₂ a year which could be incorrectly allocated depending on the Land trapping relationship extracted from core-scale experiments. Whether the CO₂ is trapped by pore-scale residual trapping or larger scale capillary heterogeneity trapping has important implications for site operators considering storage capacity and security of a site.

6 Conclusion

Capillary heterogeneity trapping has been directly observed in a heterogeneous Bentheimer core, with a local increase in non-wetting phase saturation upstream of the capillary pressure barrier. Experimentally, the proportion of capillary heterogeneity trapping is shown to be inversely proportional to brine imbibition rate. The observed initial-residual trapping relationships behaved similarly to literature relationships for homogeneous or heterogeneous cores, depending on flow rate.

The analytical model outlined can be used identify the key parameters influencing capillary heterogeneity trapping. The experimental data is observed to match with analytical model predictions. In future work, the analytical model could be used to estimate the proportion of capillary heterogeneity trapping within heterogeneous core-floods. The dimensionless trapping length is identified as the key dimensionless number to describe capillary heterogeneity trapping. In addition, by matching the analytical solution to measured initial-residual trapping characteristics, further details of the core may be obtained such as the capillary pressure of the layers or even the parameterisation of relative permeability.

Further work is needed to develop a parameterisation of the initial-residual saturation relationship that can incorporate the influence of rate on sub-core scale capillary heterogeneity trapping and can also be used in upscaling to the field scale.

Open Research

Data and software (analytical model, data analysis, data visualization, and model output) associated with this work may be obtained from the BGS data repository, item ID 187699 [Harris2025][Unrestricted access, Open Government Licence (OGL)].