Small‐Scale Capillary Heterogeneity Linked to Rapid Plume Migration During CO2 Storage

Unpredicted, rapid plume elongation has been observed at subsurface CO2 storage projects worldwide, exemplified by the Sleipner project. We show that conventionally ignored centimeter‐meter scale heterogeneity in capillary pressure characteristics can manifest as rapid field‐scale, decameter‐kilometer, plume migration. We analyze the effect in the Goldeneye field, UK, a proposed storage site with a unique combination of sample/data accessibility and generality as an archetype sandstone reservoir. We overcome previous barriers by characterizing in greater detail over larger scales—the 65 m reservoir height at cm‐m resolution—and through use of an upscaling scheme which resolves small‐scale heterogeneity impacts in field‐scale simulations. These models reveal that significant early time retardation of buoyantly rising CO2 plumes is followed by rapid migration under the caprock in the presence of anisotropic, layered heterogeneities. Lateral migration speeds can be enhanced by 200%, placing first‐order controls on fluid flow and providing a mechanistic explanation for field observations.


Introduction
CO 2 plume migration remains difficult to accurately simulate and predict at industrial-scale, subsurface storage projects. Plume elongation and early breakthrough at monitoring wells has occurred at a number of projects of various sizes, see Figure 1a. This includes Sleipner, North Sea, Norway (Chadwick et al., 2009;Williams & Chadwick, 2012); In Salah, Algeria (Ringrose et al., 2013); Frio stages I and II projects, Texas, USA (Daley et al., 2011;Hovorka et al., 2006); Cranfield, Mississippi, USA (Hosseini et al., 2013); and the Otway Stage 2 project, Otway, Victoria, Australia (Dance et al., 2019). In all of these examples, the seismically imaged CO 2 plume has extended laterally away from injection points along pathways and at rates neither predicted nor easily matched by conventional approaches to reservoir simulation. The Quest project in Alberta, Canada, provides a notable exception where lateral plume migration is significantly less than anticipated from prior modeling (Halladay et al., 2018).
A number of hypotheses have been proposed to explain the unexpected behavior. Most frequently, studies invoke the existence of a high-permeability channel or fracture, aligned with the flow direction (Cowton et (Cowton et al., 2016;Dance et al., 2019;Hosseini et al., 2013;Ringrose et al., 2009). (b) Field-scale capillary numbers associated with a highly permeable aquifer (i.e., typical of the Darcy magnitude permeability of the formations at Sleipner, Otway, Frio, and Quest) as a function of injection rate and radial distance from the injection well. In this conservative example, capillary equilibrium (N c < 1) is established after a radial distance of 47 m (Q = 1 Mt/yr), 21 m (Q = 0.5 Mt/yr), and 8.7 m (Q = 0.25 Mt/yr), respectively, over a 1 m capillary heterogeneity of magnitude ΔP c = 5 kPa. Further parameters used for the calculation are detailed in the supporting information (SI). caprock topography, this can place a leading order control on the migration pathway (Cowton et al., 2018;Gasda et al., 2013;Nilsen et al., 2017). Thermal effects and the uncertainty in the fluids' thermophysical properties (including contamination) have also been suggested as an underlying cause (Hodneland et al., 2019;Williams & Chadwick, 2017). Even accounting for the large uncertainties in subsurface system properties, precise matches to plume behavior have only been achieved using model parameters, for example, permeability, outside the range of values constrained by observations. This difficulty suggests a gap in the physics represented in the models used to simulate flow behavior. A number of studies have demonstrated an important control of small-scale heterogeneities in capillary pressure characteristics giving rise to upscaled impacts on multiphase flow (Braun et al., 2005;Cavanagh & Haszeldine, 2014;Ferrand & Celia, 1992;Gershenzon et al., 2015;Li & Benson, 2015;Park et al., 2016;Saadatpoor et al., 2010;Trevisan et al., 2017;Yeh et al., 1985). In the current context, we consider small-scale heterogeneities as lenses, laminae or bedding heterogeneities ranging from 0.1 m to several meters in thickness (Gershenzon et al., 2014;Ringrose et al., 1993), which are at least an order of magnitude smaller than the plume migration. In application to CO 2 storage, the impacts have been shown to play a significant role in inhibiting the upward migration of injected CO 2 , and in the enhancement of trapping. Past work on nonaqueous phase liquid (NAPL) contamination has shown that these heterogeneities can also result in significant enhancements in lateral migration; with NAPL, this can occur with a sinking dense contaminant plume (Braun et al., 2005). This impact on lateral migration is unexplored for CO 2 storage and is the focus of this work. Moreover, the evaluation of the fine-scale controls of multiphase flow heterogeneity have been limited by both a lack of the high spatial resolution petrophysical data from field locations required to constrain such hypotheses, and the implementation of upscaling approaches that honor the impacts of these properties at field scales. In this work we combine laboratory characterization of a case study reservoir with a rigorous multiscale upscaling approach.
We show that the impact of small, dm-m scale, heterogeneity in capillary pressure characteristics, can manifest as rapid field-scale, decameter to kilometer, plume migration. We overcome previous gaps by addressing Geophysical Research Letters 10.1029/2020GL088616 issues of both characterization and modeling. We perform a high spatial resolution laboratory characterization of the entire 65 m height of a target CO 2 storage site, the Goldeneye field in the Captain Sandstone of the UK Central North Sea (Shell U.K. Limited, 2015). This is then used in the creation of numerical models in an upscaling workflow that allows us to evaluate the kilometer-scale impacts of the fine-scale heterogeneity in a model of the Goldeneye field.

Scaling Analysis for Small-Scale Heterogeneities
First we demonstrate the general prevalence of the impacts of small-scale heterogeneities by consideration of the forces driving fluid flow. Subsurface CO 2 migration is driven by gradients in fluid pressure, capillary pressure (P c = P CO 2 − P w ) and buoyancy (Zhou et al., 1994). Conventional scaling analysis of these potentials identifies that the viscous or buoyant terms dominate at spatial scales of order 1 m and greater. Capillary pressure gradients can be ignored over field scales because they contribute little to fluid movement (Bear, 1989;Zhou et al., 1994).
This approach overlooks a key impact that small-scale heterogeneity in capillary pressure characteristics can have on large-scale plume migration. The establishment of equilibrium in capillary pressure over small length scales in heterogeneous rocks has a significant impact on the local fluid mobility. As has been identified previously, trapping will be significantly enhanced with fluid flow across layers (Krevor et al., 2011;Meckel et al., 2015;Saadatpoor et al., 2010). However, average relative permeability can also be dramatically enhanced in flow parallel to layers across which capillary equilibrium has been established (see Virnovsky et al., 2004, for numerical simulations, andKrevor, 2015, for experimental observations).
We illustrate the prevalence of the state of local capillary equilibrium with a simple scaling argument. A dimensionless capillary number, N c , represents the ratio of viscous to capillary forces. We use the number proposed by (Virnovsky et al., 2004), N c = HΔP LΔP c , but the discussion applies using equivalent numbers proposed by others, for example, (Zhou et al., 1994). There are two key spatial scales. The flow length L is the path length along which fluid is driven by buoyancy or gradients in pressure, ΔP. Second, the spatial scale H is the length-scale characteristic of a heterogeneity, for example, the thickness of a bedform layer in sandstone. This length scale is limited to the distances, approximately meters or less, over which fluid distribution by capillary forces is approximately linear with respect to the gradient in capillary pressure. Capillary equilibrium is established when gradients in capillary pressure, ΔP c , equalize across heterogeneities. The flow path length, L, at which the capillary number N c < 1 indicates the distance of fluid travel, for example, away from an injection point, beyond which capillary equilibrium across heterogeneities has been established. This is the point at which enhancements to fluid mobility and trapping due to heterogeneity will prevail.
Approximations of these length scales for existing injection sites shows that capillary equilibrium is generally established within several meters of a CO 2 injection well. This is due largely to the radial nature of the injection, meaning the driving pressure decays logarithmically away from the source; even in the strongest injection settings, capillary heterogeneities will eventually prevail. Figure 1 provides the estimated equilibrium length scale (in this case the radial distance, r, is used, that is, L = r) for varying injection rates for an archetype, highly permeable aquifer, analogous to the Darcy magnitude sandstones at Sleipner, Frio, Quest, and Otway, with calculation details in the supporting information (SI). This example is conservative, with lower permeability, stronger heterogeneity contrast cases (e.g., In Salah and Cranfield) likely to equilibrate within shorter distances.
Thus, it is a general feature of subsurface CO 2 migration over large scales that it is characterized by flow where the fluid distribution is in capillary equilibrium across heterogeneities up to meters in dimension. Fluid mobility (e.g., Kk r / ) will therefore be largely controlled by the nature of small-scale heterogeneities in multiphase flow properties; flow along layers will have enhanced mobility, while flow against layers will have reduced mobility and increased trapping. In the following sections we evaluate the quantitative impact these small heterogeneities have on plume migration away from injection points.

The Impact of Small-Scale Heterogeneities on Lateral Plume Migration
We now evaluate the impacts of capillary equilibrium across small-scale heterogeneities on lateral plume migration in high spatial resolution simulations in which we explicitly resolve the small-scale features. This is done to observe the impacts without incurring numerical errors associated with upscaling, and also as  (Shell data, Shell U.K. Limited, 2015) against rock core plug sample data from nitrogen porosimetry (Imperial College London (ICL) regular core analysis laboratory (RCAL) data). The solid gray line represents the porosity depth trend = −0.0003TVD + 1.01, with dashed lines showing the standard deviation measured from combined Shell and ICL data. (b) Porosity-permeability relationship, the solid black line represents ln(K) = 52.544 − 6.4656. (c) Mercury intrusion capillary pressures for 13 plugs comprising the composites, scaled using the measured porosity and permeability (Leverett, 1941). (d) Semivariogram plot of porosity at increasing lag distance (depth). (e) Experimental relative permeabilities for composite cores through the reservoir. Core depths: Composite core 1-2552.4-2555.4m, Reynolds et al. (2018Reynolds et al. ( )-2562Reynolds et al. ( .7-2568, Composite core 2-2603.6-2604.9m. a first step in developing an upscaling approach. We use a geological setting that was the focus of significant appraisal and development activity for an industrial-scale CO 2 storage project-the Captain Sandstone Formation in the Goldeneye field, North Sea, UK (Tucker & Tinios, 2017). The primary injection unit, the Captain D, is a poorly consolidated, medium-grained massive sandstone which fines upward with typical permeabilities in the range 0.7-1.5 Darcy. It is an archetype unit in the North Sea and draws parallels with other active CO 2 storage sites such as Sleipner and Snøhvit on the Norwegian Continental shelf (Eiken et al., 2011).
We characterize heterogeneity in the reservoir through laboratory analysis of 46 rock samples, sampled from an exploration well (14/29a-3) comprising the entire 65 m vertical interval of the targeted injection site (summarized in Figure 2, experimental methods are detailed in the SI). The porosity and permeability are correlated and show a decreasing trend with depth (Figures 2a and 2b). The capillary pressure-saturation characteristics were scaled using the Leverett-J-scaling (Leverett, 1941), collapsing onto the same intrinsic function describing the interval. Similarly, we find little variation in the viscous limit, intrinsic relative permeability function through the formation (Figure 2e). Heterogeneity in the formation appears in the form of dispersed, low-permeability mud clasts (Shell U.K. Limited, 2015), with small layers and dish-and-pillar structures of vertical correlation length 1-4 m, visible in the core, well-log and semivariogram data, see  (b) Relative migration speed, V r , at low flow potential as a function of correlation length ratio (anisotropic cases). V r is the inverse of the relative breakthrough time-specific realization breakthrough times, t i , are scaled against the uncorrelated average for the corresponding heterogeneous or homogeneous P e case, t 0 , in which r x = r y = 0.1 m (bottom plots in (a)). Each data point represents the mean from five different geostatistical realizations, with error bars showing the range. Multiple points at a given r x /r y have different absolute r x and r y , but maintain the aspect ratio. The solid, dash-dotted and dotted lines highlight log curve fits through U d = 0.002, 0.2, and 20 m·day −1 data, respectively. Only the U d = 0.002 m·day −1 data points are shown explicitly for clarity. (c) Plume breakthrough times for heterogeneous P e cases as a function of capillary number N c . 80 • C (Shell U.K. Limited, 2015) and average hydrostatic pressure of 25.26 MPa; fluid properties, geological and simulation setup are found in the SI. We vary the correlation length of heterogeneities to generate layers of varying lateral length. Vertical correlation lengths, r y , are varied between 0.1 and 5 m, and horizontal correlation lengths, r x , from 0.1-50 m, consistent with the petrophysical evidence. Porosity, permeability and capillary pressure heterogeneity are fully correlated through the petrophysical relationships and Leverett-J scaling detailed in Figures 2a-2c, and further in the SI. The spatially homogeneous, intrinsic relative permeabilities from Figure 2e are prescribed everywhere in the domain. We simulate the isothermal, incompressible lateral migration of CO 2 and 1 Molal NaCl brine. The domain is initially brine saturated at hydrostatic pressure, with CO 2 subsequently injected at a constant rate through the left-hand boundary.
We observe up to a 200% increase in plume migration rate as a result of the layered, anisotropic heterogeneities in the capillary pressure characteristic (Heterogeneous P e in Figure 3b). The results are scaled against a base case representing the mean of the isotropic, uncorrelated cases (Figure 3a, bottom heterogeneous). We see similar qualitative enhancements when the migration speeds are compared directly to the corresponding homogeneous cases (see SI Figure S16; migration is enhanced by up to 100%). Figure 3 shows that the CO 2 plumes in anisotropic cases are elongated in the domain compared to isotropic cases and homogeneous P e cases. The plumes migrate laterally much further for the same injected volume. The CO 2 has preferentially saturated low entry pressure regions (visible in the entry pressure maps) creating high-mobility regions through which the CO 2 can rapidly migrate. This applies generally to simulation results under the range of heterogeneity realizations and injection rates shown in Figures 3b and 3c. They cover a wide range of flow regimes defined by the capillary number, N c , created by varying the lateral injection velocity U d from 0.002-20 m·day −1 , changing the force balance from capillary to viscous dominated. As the degree of anisotropy in the capillary pressure heterogeneity increases, the migration speed of the plume increases significantly (Figure 3b, isotropic cases are shown in the SI).
The impact is also almost entirely a result of heterogeneity in capillary pressure characteristics alone, and not permeability and porosity. The homogeneous results in Figure 3a and the black points in Figure 3b are results from simulations with domains heterogeneous in porosity and permeability, but homogeneous in capillary pressure characteristics. There is relatively little impact on the plume shape and breakthrough time when varying the correlation length of the permeability and porosity.
The lateral migration speed is generally faster over the entire range of capillary numbers when the dimensionless breakthrough time from the anisotropic cases are compared to the isotropic cases in Figure 3c. However, they are generally much faster at low to moderate capillary numbers, with the impact of heterogeneities becoming less pronounced at high capillary numbers. At larger correlation lengths and length ratios, the variance in breakthrough time increases, indicating the increase in flow path tortuosity and system representativeness. The error bar ranges reported in Figure 3b indicate this variance, showing that the same general trends hold at larger sizes. There are some outliers in the anisotropic data, with slow breakthrough times; these occur in cases where both correlation lengths approach the respective system size, which in some instances can retard the flow.
As N c increases in Figure 3c, the breakthrough time correspondingly increases, until it reaches a plateau at N c ≈ 10 2 . The increase in breakthrough time occurs due to the greater vertical sweep of the plume as the flow rate and N c increase, and the plume behavior in the plateau regime is usually referred to as viscous dominated. The viscous limit at N c ≈ 10 2 is similar to both previous experimental and modeling studies at a range of scales (Jackson et al., 2018;Virnovsky et al., 2004). The N c reported in Figure 3c also highlights the increase in pressure buildup in the isotropic systems compared to anisotropic systems, for a given flow rate. For a fixed pressure buildup, the isotropic heterogeneity reduces the overall system injectivity. The elevated pressure plume is also a concern in the far field, as it diffuses and extends beyond the physical CO 2 footprint, potentially causing brine leakage into shallow formations among other effects (Cihan et al., 2013).

Field-Scale Implications and Discussion
Conventionally, the field-scale impact of capillary pressure has been evaluated without consideration of small-scale heterogeneities. In homogeneous reservoirs capillarity leads to an increase in the sweep of the plume and retardation of its lateral migration (Becker et al., 2017;Golding et al., 2013). We have shown that the presence of structured heterogeneities can lead to an outcome which has qualitatively the opposite impact. Capillary heterogeneity provides a mechanism for rapidly elongating and accelerating the plume migration. We now demonstrate the field-scale implications of the above results by considering a full-field model of the Captain sandstone formation in which rigorous upscaling techniques have been applied to honor the impacts observed at high resolution.
We use a novel upscaling approach (detailed and validated in the SI) to demonstrate the combined impacts of layered capillary heterogeneities in a field-scale setting, during both vertical and lateral migration regimes. The upscaling method allows the incorporation of fine-scale capillary heterogeneity impacts in models with much larger cell sizes, through generation of equivalent, macroscopic flow functions. This permits analysis of large systems that would be insurmountable with fine-scale grids.
We use a 2 km east-west transect through the center of well 14/29a-3 to generate a 2-D model of the Captain D sandstone in the Goldeneye field. The model has capillary pressure heterogeneities in line with previous petrophysical analysis, with correlations r x = 50 m, r y = 2 m. Carbon dioxide is injected with outlet velocities of 0.01 m day −1 into a central well (larger injection rates are shown in the SI), which is perforated through the bottom one third of the domain, shown in Figure 4a.
Snapshots of the plume location at different times are shown in Figures 4b-4d, illustrating the key impacts of small-scale capillary heterogeneities in field-scale CO 2 flow. The top plots show the case with capillary heterogeneity, the bottom plots show the homogeneous case (still with permeability and porosity heterogeneity to isolate the impact of capillary heterogeneity). Alongside this, the distance from the top of the plume to the caprock (R 1 ) and the lateral migration distance of the plume along the caprock (R 2 ) are shown through time in Figure 4e.
The initial buoyant rise of CO 2 is significantly retarded by the anisotropic capillary heterogeneities, which act as vertical flow baffles. The flow is spread underneath these baffles, pooling until the column height underneath is great enough to overcome the high capillary entry pressure of the local heterogeneity. The migration time to reach the cap rock is 0.035 pore volumes for the heterogeneous case, but only 0.005 pore volumes for the homogeneous case, a sevenfold decrease.
The impact of heterogeneities on lateral migration is observed after the initial buoyant rise. The CO 2 sweeps laterally underneath the cap rock, still largely driven by the buoyancy of the vertical CO 2 column underneath, shown in Figures 4c and 4d. In the heterogeneous case, the plume extends more thinly underneath the caprock, migrating in regions of low entry pressure. Indeed, given the initial retardation of the vertical migration, the plume migrates much quicker laterally in the heterogeneous case compared to the homogeneous case, as can be seen in Figure 4e. It migrates in high-speed "bursts" when the flow path aligns with the anisotropic heterogeneities, ultimately catching the homogeneous plume. The results here show the combined impact of layered capillary heterogeneities in the field when flow migration is aligned predominantly perpendicular and parallel to layers, respectively. The inhibiting impact of layering on upward migration and trapping has been studied previously (Li & Benson, 2015). This effect is dependent on the extent to which the layers cannot be bypassed by the fluid. Layers with shorter correlation lengths in one lateral dimension will provide a bypass route for the flow, partly offsetting the inhibition effect.
The case herein represents a relatively conservative scenario for the enhancement of lateral migration speed due to the topography of the Goldeneye field. The steep anticline causes CO 2 to pool and migrate both laterally and downward against the heterogeneities, mitigating the enhanced rate of plume migration. In cases of a flat topography (i.e., in the examples in section 3), an updip migration, or a shallow anticline (i.e., at Sleipner, Furre et al., 2017), the enhancement would be significantly more prominent. Although the final lateral extent appears similar in Figure 4, the lateral spreading is much faster in the heterogeneous case ( Figure 4e).
Small-scale, capillary heterogeneities could have contributed significantly to the rapid plume migration, elongation, and early monitoring well breakthrough that has been observed in several field sites, that is, Sleipner, North Sea, Norway (Chadwick et al., 2009;Williams & Chadwick, 2012); In Salah, Algeria (Ringrose et al., 2013); Frio stages I and II projects, Texas, USA (Daley et al., 2011;Hovorka et al., 2006); Cranfield, Mississippi, USA (Hosseini et al., 2013), and the Otway Stage 2 project, Otway, Victoria, Australia (Dance et al., 2019). Small-scale capillary features (i.e., 1-10 s of meters) are generally not included in the field-scale simulation of these sites, largely due to model and/or geological resolution. The leading order impacts of small heterogeneities suggests that these problems should be revisited making use of characterization, upscaling, and modeling techniques that honor the impacts these heterogeneities have on field-scale flow.
The characterization and upscaling methodology herein provides a method for tractably simulating large-scale flows at low capillary number, including the impacts of small-scale heterogeneities. This is particularly important in uncertainty analysis, allowing different geological realizations to be run quickly to see potential impacts of different interpretations, and when new information comes through from field operations (Dance et al., 2019). The ability of numerical models to conform to observed plume migration is vital for demonstrating compliance with regulators, that is, within a specified area of review (Court et al., 2012;