The impact of capillary heterogeneity on CO₂ plume migration at the Endurance CO₂ storage site in the UK Open Access

• Upscaling the impact of capillary heterogeneity on the multiphase flow properties is important for modelling subsurface CO₂ flow.

• Capillary heterogeneity significantly enhances preferential flow pathway formation, increases reservoir sweep and reduces vertical migration speed near CO₂ injectors.

• Simulations in 3D exhibit qualitatively distinct results from those in 2D due to the increased tortuosity in flow pathways.

Numerical simulation is used to support development planning and uncertainty analysis in geological CO₂ storage projects. Models provide insights into volumes stored, optimum well locations and bounding risks associated with geological uncertainty. However, several field scale projects have reported plume dynamics that do not match numerical predictions, most notably at Sleipner (Chadwick and Noy 2010). This suggests that some important flow behaviours may be missing from field scale flow simulations.

Heterogeneity in the multiphase flow properties has been identified as an important control on CO₂ plume migration, that is typically overlooked in conventional modelling workflows (Chang and Yortsos 1992; Krause 2012; Jackson et al. 2018). Permeability and porosity heterogeneity can be observed at a range of scales (⁠$μ$m - 100 s m). These heterogeneities are also associated with spatially varying capillary pressure, termed capillary heterogeneity. Studies have found that capillary heterogeneity is one of the dominant fluid distribution mechanisms over centimetre to metre length scales (Chang and Yortsos 1992; Ringrose et al. 1993; Chaouche et al. 1994; Kjonsvik et al. 1994; Huang et al. 1995; Krause 2012; Jackson et al. 2018). Further studies have demonstrated that it may also influence field scale behaviours (Kjonsvik et al. 1994; Jackson and Krevor 2020; Boon and Benson 2021; Ni et al. 2021) potentially modifying the footprint, vertical and horizontal migration distance and migration speed of the CO₂ plume (Al-Khdheeawi et al. 2017; Onoja and Shariatipour 2019; Jackson and Krevor 2020; Shao et al. 2022). The trapping of CO₂ and the storage capacity are also strongly influenced by the multiphase heterogeneity (Hovorka et al. 2004; Ambrose et al. 2008; Krevor et al. 2015; Gershenzon et al. 2017; Ni and Meckel 2022; Bump et al. 2023).

Various works, including Jackson and Krevor (2020); Ni et al. (2021); Shao et al. (2022) modelled the field scale impacts of sub-metre scale capillary heterogeneity, but its impact on the performance of a 3D, multi-well, industrial scale project has yet to be investigated. The reservoir models in Jackson and Krevor (2020) and Ni et al. (2021) were limited to only two dimensions, thus they were not able to consider the changing importance of capillary, viscous and gravity effects with distance from the injection well. Shao et al. (2022) simulated 3D plume migration at a pilot site, however this had low injection volumes, a single well in the centre of the model and a limited domain size.

In this work we study the impact of heterogeneity on plume migration at an industrial scale geological carbon storage site with multi-point CO₂ injection over a 25 year time scale. For this purpose, we focus our investigation on the proposed Endurance CCS site, which is characterized by strong lateral and vertical rock heterogeneity (Gibson-Poole et al. 2024). The site is a four-way dip closure, which enables us to investigate the interplay between buoyancy and heterogeneity in more detail. We upscale the impact of metre scale capillary heterogeneity on the >100 m scale multiphase flow properties, the capillary pressure and the relative permeability. We simulate kilometre scale CO₂ plume migration under the influence of capillary heterogeneity and explore the impact of well location on the resultant plume migration. Figure 1 summarizes the flow behaviour discovered in this work, which is attributed to the presence of capillary heterogeneity.

We construct a model representative of the central part of the Endurance anticline, the storage reservoir for Phase 1 of the Northern Endurance Partnership (BEIS 2021c). The site is characterized by strong heterogeneity in all three dimensions with a distinct, four-way dip closure. The injection wells are located downdip of the structural high, thus the CO₂ will tend to migrate up-dip away from the injection wells. The flow will be influenced by gravity as well as capillary heterogeneity and viscous effects. The model has overall dimensions (x, y, z) 6000, 4000 and 250 m, and encloses the approximate location of four out of five CO₂ injectors planned for phase 1 of the project, in the central portion of the Endurance anticline (Fig. 2).

A fine grid model of the central part of the Endurance anticline was used to determine the upscaled petrophysical properties, with grid block dimensions (x, y, z) of 10, 10 and 1 m. The petrophysical properties were derived from the well log data, coreflood measurements and geological constraints provided in BEIS (2021c). Their distribution was modelled using the geostatistical approach described by Jackson and Krevor (2020), with inputs taken from characterization studies performed on the Endurance site (BEIS 2021c). We opted for a geostatistical approach because it allows us to clearly identify the role of salient features including correlation lengths and degree of heterogeneity on plume migration. It is also consistent with the geological modelling approach described in BEIS (2021c).

A heterogeneous, spatially correlated porosity field was generated using ellipsoidal averaging (Durlofsky et al. 1997), following the approach applied to a 2D model by Jackson and Krevor (2020). The standard deviation was obtained by fitting an extreme value distribution to the porosity data of well 42/25d-3 (Fig. 3c). The correlation lengths in our model were set to (r_x, r_y, r_z) 1000, 760 and 4 m. The horizontal correlation lengths, r_x and r_y, were based on the data provided in BEIS (2021c). The vertical correlation length, r_z, was obtained by generating a semi-variogram using the porosity data from well 42/25d-3, as shown in Figure 3d. The initial sill distance of 4 m, plotted as a purple line in Figure 3d, indicates the correlation length of the vertical heterogeneity. The variogram is also characterized by periodicity at larger lag distances, which is indicative of depositional periodicity (Gringarten and Deutsch 2001; Jackson and Krevor 2020).

The correlated field was re-scaled to match the mean and standard deviation from the extreme value model of the porosity well data, and the depth transform from Figure 3a was applied. This resulted in a correlated field that matches the statistical information obtained from characterization and data collection studies conducted as part of the Northern Endurance Partnership project.

After the correlated field was generated, the topography of the four-way dip closure, characteristic of the Endurance site, was created. The top surface of the initially cuboid fine grid model was modified with a depth transform that follows the topography of the top Bunter sandstone presented on contour maps in White Rose (2016). The resultant top structure of our model is presented on a topography contour map in Figure 4.

The fine-scale grid block permeability was calculated using the porosity-permeability trend from Figure 3b. We derived the fine-scale entry pressure distribution using J-function scaling of the capillary pressure characteristic curves collected using mercury intrusion porosimetry on samples from the Endurance site.

The fine-scale model has 60 million grid blocks. In order to tractably simulate multiple injection scenarios, we upscaled the fine-scale grid. A capillary-limit upscaling approach was used to determine the multiphase flow properties.

We based our approach on the work by Jackson and Krevor (2020), who developed a steady-state, capillary-limit upscaling scheme using macroscopic invasion percolation on a fixed, uniform Cartesian mesh for 2D domains. Their method is grounded in the theoretical developments of macroscopic invasion percolation and steady-state upscaling, primarily from Wolff et al. (2013) and Yang et al. (2013), but also from Pickup and Stephen (2000), Virnovsky et al. (2004), Braun et al. (2005), Odsater et al. (2015), Hilden and Berg (2016) and Rabinovich et al. (2016). We extended this workflow to the third dimension and applied it to the fine-scale reservoir model of the Endurance site.

The Endurance site is characterized by heterogeneity in all three dimensions indicated by the correlation lengths implemented in the geological model. Endurance also has a distinct, four-way dip closure, which will likely interact with the CO₂ as it migrates through the pore space. Therefore, extending this upscaling scheme to 3D allows us to understand the leading controls on the CO₂ flow behaviour. The approach and verification is described in the supplementary information (SI). The upscaled model has grid block dimensions (x, y, z) of 200, 200 and 3 m. The upscaled heterogeneous porosity, permeability, and entry pressure are displayed in Figure 5.

We used a commercial, black oil simulator (Schlumberger's E100™) to model CO₂ injection in our upscaled reservoir model. We assumed the reservoir was isothermal as the focus of this work was on capillary heterogeneity upscaling and its impacts on plume dynamics. Temperature variations may additionally impact CO₂ migration as shown by Al-Khdheeawi et al. (2018) and Abdoulghafour et al. (2020). In the case of the Endurance site, observations have been reported from the 42/25d-3 appraisal well of a 7°C variation across the 200 m vertical height of the reservoir (White Rose 2016). These will impact the plume primarily through an impact on the viscosity of the CO₂. However, the variation is modest and small differences in the distribution of the plume will not have a qualitative impact on the findings.

The reservoir temperature was set to 56°C (White Rose 2016) and the initial reservoir pressure was determined assuming hydrostatic equilibrium, with the depth of the top of the anticline set to 1020 m. This results in an average initial reservoir pressure of 146.4 bar. The fluid properties under these conditions are provided in the SI.

Another key simplification is that we did not simulate CO₂ dissolution. Studies including Lindeberg and Begmo (2003), Ampomah et al. (2016), Gilmore et al. (2020) and Yuan et al. (2020) suggest that dissolution may be an important mechanism to consider in CO₂ injection simulations, particularly in the presence of heterogeneity. However, the upscaling of coupled physical processes is complex, particularly for CO₂ dissolution, which may require consideration of multiple types of solute transport, e.g. diffusion, dispersion, advective mixing. The focus of this work is thus to illustrate the impacts of heterogeneities on the free phase flow of CO₂.

We also implemented grid refinement around the wells to limit numerical dispersion. The details of this are discussed further in the SI. In the analysis of the simulation results, we set a CO₂ saturation threshold of S_CO2 = 5 × 10⁻³ to remove noise and numerical uncertainties from the results. This threshold is able to detect 95% of the CO₂ plume at early injection times (days-weeks), which increases to 98% for later times (months-years).

In this section, we describe the results from simulating CO₂ injection in the full 3D model at a constant rate of 3131 m³d⁻¹, corresponding to 0.8 Mt yr⁻¹, from well CI2 (Fig. 2). This rate matches the average injection rate proposed in the field development plan (BEIS 2021a). Rather than injecting from all four planned wells, we initially used a single well to remove the complexities of multi-plume interactions. We increase the number of wells in the following subsection. Constant pressure boundary conditions, which are widely implemented in CO₂ injection simulation studies (Doughty 2010; Zhou and Birkholzer 2011; Shao et al. 2022), were implemented on all four sides of the domain. This allows fluids to leave the system, which is a suitable choice considering our model domain is smaller than the entire Endurance structure.

We generated six models using a combination of data from the Endurance site to differentiate the impact of capillary pressure variations and relative permeability anisotropy on CO₂ flow behaviour. The six models can be categorized into two groups: one that incorporates the heterogeneous and anisotropic relative permeabilities that are the result of upscaling (referred to as ‘CL’ for capillary limit) and one that assumes constant relative permeability set to the fine-scale viscous-limit relative permeability (referred to as ‘VL’ for viscous limit). For each category, we created three models: one with the heterogeneous, upscaled capillary pressure (‘Pc’), one with constant capillary pressure set to the entry pressure of the fine-scale P_c–S_w curve (‘Pe’) and a third with capillary pressure set to zero (‘Zero’). Thus the six cases from most to least heterogeneous are:

• CL Pc: heterogeneous, upscaled relative permeability and capillary pressure

• CL Pe: heterogeneous, upscaled relative permeability and constant capillary pressure

• CL Zero: heterogeneous, upscaled relative permeability and zero capillary pressure

• VL Pc: homogeneous relative permeability and heterogeneous, upscaled capillary pressure

• VL Pe: homogeneous relative permeability and constant capillary pressure

• VL Zero: homogeneous relative permeability and zero capillary pressure

Figure 6 shows the upscaled block relative permeabilities and capillary pressures from the y–z plane at x = 20 in the 3D model. We observe an increase in the upscaled x and y gas relative permeabilities compared with the intrinsic curves (Fig. 6b, c), which is a result of the large correlation lengths in the x and y dimensions. This may enhance gas flow in these dimensions. In comparison to this, the majority of blocks exhibit a decrease in the vertical gas relative permeability (Fig. 6a), which can also be attributed to the heterogeneity structure in the model. For all dimensions, there is a wide range of gas relative permeability curves. This is a result of the adverse viscosity ratio of CO₂ displacing brine combined with the heterogeneity correlation length being much greater in the x and y directions than z direction. We also observe a higher trapped gas saturation for the z direction saturation functions, than x or y. This a consequence of the vertical correlation length being similar to the coarse grid block size in the z direction resulting in increased capillary trapping of the CO₂ at layer boundaries.

To analyse the variation in plume migration speed between the six simulations, we extract the maximum spread of CO₂ in the x, y and z dimensions as a function of time. This is defined as the difference between the maximum and minimum coordinates in a single dimension invaded by CO₂ (S_CO2 >5 × 10⁻³). The steeper the gradient, the faster the plume is travelling in that dimension. This is presented in Figure 7.

From the onset of injection, the plume migration speed in the x-direction is higher for the models with heterogeneous, upscaled relative permeabilities (plotted as purple lines in Fig. 7a). The difference in x lateral migration between the CL and VL models increases with time; the plume has spread an additional 600 m after 25 years of injection, which represents 10% of the total domain size in the x dimension. Thus, capillary heterogeneity enhances plume spreading in the x-direction.

The trend in the migration speed in the y-axis is slightly different (Fig. 7b). Initially, the CL models are characterized by higher speeds. However, the plumes in the VL simulations soon overtake and migrate faster between 1.7–6.8 years. This is followed by a region where the plumes in the CL models spread further. After 25 years, the CL models have spread an additional 200 m in the y-axis. By comparing the observations from Figure 7a and b, it emerges that all plumes generally spread further along the x-axis than the y-axis. This is partly influenced by the topography implemented in the model. The topographic high is located at grid block location x = 19, y = 12. Thus, buoyancy will result in the plume migrating more in the x than y direction. The correlation length of the heterogeneity also affects plume shape. In our model, r_x is 30% larger than ry. Therefore, the layers are more continuous in the x-axis, increasing the tendency of the plume to migrate in that direction.

The enhancement of the lateral migrated distance and lateral migration speed due to capillary heterogeneity is not observed in the 2D simulations presented in the SI. Instead, the plume in the 2D VL models travelled nearly 1 km further laterally over the 25-year period. These results imply that the impact of capillary heterogeneity on the lateral migration of the plume is qualitatively distinct in the 2D models. One reason for this may be the difference in the heterogeneity structure imposed in the 2D v. 3D domains. In 3D, the third dimension opens the possibility for many alternative flow pathways, which can help the plume avoid capillary barriers. In 2D this effect is lost and the plume may be held up behind capillary barriers (low k_r/high P_e barriers in the CL Pc simulation, low k_r barriers in the CL Zero simulation), which results in a less efficient migration. This shows a major qualitative difference to the 3D models. The reader is referred to the SI for a more in-depth discussion.

Up to 5.4 years, there is a higher vertical migration speed in the VL models, reaching a maximum difference of 80 m from the injection point to the highest z location. Ultimately, the faster vertical migration means the plume reaches the caprock after only 74 days in the VL cases, whereas it takes nearly 400 days longer in the CL runs to do so. This is consistent with observations from 2D simulations discussed in the SI. After 5.7 years the behaviour changes and there is a larger vertical spread in the CL models, in other words the plume in regions beyond the injection point spreads over a larger vertical distance in the CL models. Reasons for this are investigated in subsequent sections.

Interestingly, the impact of capillary pressure variations on migration speed in any dimension is negligible for both the CL and VL runs. This indicates that with this combination of coarse grid block dimensions and heterogeneity correlation length, the relative permeability functions and their anisotropy, are the main control on migration speed. It is important to remember that the fine-scale capillary pressure affects the shape of the upscaled relative permeability curves, and thus its impact is represented with these flow functions. The greater importance of the upscaled relative permeability curves relative to the upscaled capillary pressure curves is consistent with scaling theory, which shows that viscous forces dominate over capillary forces at greater length scales (Jackson and Krevor 2020). For most reservoir systems, these length scales are achieved at 10 m or greater.

Additional comparisons between 2D and 3D simulations presented in the SI indicates that upscaled capillary pressure has a more measurable impact in 2D (although upscaled relative permeability is still the dominant factor). This suggests the main role of the upscaled capillary pressure function in the simulations is through the capillary entry pressure, which excludes CO₂ from grid blocks where that has not been exceeded. In 2D, the CO₂ can be trapped upstream of these locations more readily than in 3D, where there are more pathways around the obstructions.

The arrival time of the CO₂ at the top structure is 8.5 years for the VL Zero simulation, while the plume in the CL Pc model arrives roughly 2.5 years later. The slower vertical migration in the CL simulations at early times influences the lateral distance travelled by the plume at later times. These results imply that small-scale capillary heterogeneity has a stronger upscaled impact on total plume migration speed.

We define plume footprint area as the area of the CO₂ plume when it is projected onto a horizontal plane. This captures plume spreading in any vertical layer of the reservoir. It would be positively correlated with the amount of CO₂ that is immobilized due to dissolution and capillary trapping, if these were modelled (Sundal et al. 2015). The plume footprint area over the 25 year injection period is plotted in Figure 8 for the six simulations. The differences between the VL and CL runs can be categorized into four distinct regions, which are annotated accordingly. We will discuss each of these regions separately in the following.

At early times (region A in Fig. 8), the footprint area of all the cases is very similar although the CL cases have a slightly higher footprint area compared with the VL cases. The vertical distribution of the CO₂ is however very different between the CL and VL cases. In the VL case, the CO₂ primarily rises upwards and accumulates under the caprock (Fig. 9). In contrast, the vertical heterogeneity hinders the buoyant rise of CO₂ in the CL run, resulting in the CO₂ moving laterally through various layers near the perforation interval (Fig. 10).

The relative percentage difference in footprint area between the CL and VL runs decreases from 400% to 20% within the first year, until the curves crossover at 1.7 years. The VL cases now have a larger footprint area (region B in Fig. 8), which is a result of the plume having risen rapidly due to buoyancy and then spread underneath the seal. The CO₂ in the VL runs reaches the caprock 400 days earlier and thus begins migrating laterally underneath the caprock confined to the uppermost layers. The larger footprint area in the VL models during this time period is attributed to a larger plume spread in the y direction, as shown in Figure 11.

At 4.4 years, the footprint areas converge for all six cases, as indicated by annotation C in Figure 8. The results displayed in Figure 11 show that both plumes are spreading out laterally in this period. Interestingly, the plume in the VL simulation is more circular, extending a similar distance in x and y, whereas the plume in the CL run trends towards the x-axis. This is in agreement with the results displayed in Figure 7.

Beyond roughly 6.9 years, the CL models have a larger footprint area, particularly due to the increased spread in the x-direction. In this time period, we also observe a slight deviation in footprint areas between the three capillary pressure cases, which suggests that capillary pressure variations become more important at later stages. This is in agreement with findings from 2D simulations discussed in the SI.

The plume swept volume is defined as the volume of all grid blocks invaded by CO₂ (S_CO2,block >5 × 10⁻³). The swept volume provides an indication of the efficiency of pore space use, which is directly linked to the CO₂ storage efficiency of a reservoir. It also gives further insight into more subtle differences in the plume shapes, particularly the formation of thin preferential flow pathways that may not be noticeable when considering footprint area.

The plume swept volume for the six simulations is shown in Figure 12. From the onset of injection, all three CL models are characterized by a larger swept volume, which increases significantly with time. The largest difference occurs at 1.6 years (130% relative difference). These results demonstrate that reservoir sweep is significantly enhanced in the presence of capillary heterogeneity. After approximately 2 years, we also observe a deviation between the three upscaled capillary pressure cases (Pc, Pe and Zero) of the CL runs, which enhances with time. This indicates it is important to include the full upscaled capillary pressure curves in the model, not just the upscaled relative permeability curves and the capillary entry pressure.

To investigate the impact of capillary pressure variations in more detail, we extract the CO₂ saturation maps for the CL Pe and CL Pc simulations (Fig. 13). The saturation distribution for the CL Zero run is not shown as this behaves similarly to the CL Pe model. As indicated by the plume swept volume, a higher number of preferential flow pathways form in the model with upscaled capillary pressure, particularly in proximity to the injector. These results indicate that it is the variations in capillary pressure with saturation that affect the CO₂ flow pattern, not just the capillary entry pressure.

The cross-sectional saturation maps also highlight the formation of a laterally extensive preferential flow pathway that is orientated parallel to the x-axis and is present in both the CL Pc and CL Pe simulations (see annotation in Fig. 13). This preferential flow pathway does not form in the VL simulations, which suggests that it is caused by capillary heterogeneity. The upscaled x relative permeability curves of grid blocks in proximity to this preferential flow pathway are plotted in Figure 6b. As shown, the capillary heterogeneity has resulted in raised relative permeabilities, which are favourable flow pathways for the CO₂. Preferential flow paths are further enhanced by the adverse mobility ratio between the CO₂ and brine. These findings are in agreement with previous works including Chang et al. (1994), Chasset et al. (2011); Hesse et al. (2006) and Shao et al. (2022). However, our work demonstrates that it is both, the anisotropic relative permeability and the capillary pressure variations, that cause the prefernetial flow paths for the CO₂. The SI presents further analysis, in which we calculate the difference in plume thickness between the VL and CL simulations at each (x, y) location at different times. The results show that the CL Pc simulation has an increased sweep efficiency near the CO₂ injector.

We increase the level of complexity of the model to more closely replicate the Endurance field development plan. Rather than injecting from one well, we now inject from all four wells that are enclosed in our model domain (see Fig. 4). Each well injected at 3131 m³d⁻¹, which also matches the field development plan. We focus our analysis on the two end-member models: CL Pc and VL Zero.

The trends in plume area and swept volume are similar to those seen in the single well simulations (Fig. 14) suggesting that the single well results are not just due to a particular realization of heterogeneity near that particular well. The difference between the VL Zero and CL Pc footprint areas is on average 7% larger than in the single well model. This is because the plume is sampling more of the variation in the reservoir properties and caprock topography. Through this effect, the impact of capillary heterogeneity on lateral plume spreading is enhanced with four injectors. After 25 years of injection, the relative difference in swept volume between the CL Pc and VL Zero model is 38%, which is approximately 10% smaller than the difference exhibited in the single-injector scenario. This is attributed to the increase in total CO₂ volume injected into the domain. Although preferential flow pathway formation remains enhanced in the CL model with four injectors, these begin to merge over time as we introduce more gas and fill the reservoir. Ultimately, as progressively more of the reservoir is filled, the difference in the swept volume between the two simulations begins to decrease.

The effect of topography on plume behaviour is seen in the saturation distribution maps (Fig. 15). The plume at CI1 spreads out the most, whereas the plumes at CI5 and CI2 spread out the least. As more CO₂ is injected, the four plumes begin to merge. The time at which this occurs is determined by the initial vertical plume migration speed. The plumes in the VL Zero models begin to migrate laterally sooner and so they merge sooner than in the CL Pc simulations.

Interestingly, the central portion of the uppermost layers is swept last (Fig. 15c, d). This may be a result of the grid block orientation relative to the topography, which is discussed in more detail in the SI. After 7.5 years, the leading edges of the plume in the CL Pc model are less uniform than in the VL Zero case. Preferential flow pathways form in all directions from the injectors, which result in an enhanced plume spread. Preferential flow pathway formation is more closely examined in the SI using sideview figures.

These observations are consistent with those from the single-injector scenario discussed in preceding sections: capillary heterogeneity improves storage efficiency. The findings also demonstrate that the enhancement of near-injector preferential flow path attributed to capillary heterogeneity is independent of the CO₂ injector location in these reservoir models.

To further analyse the variation in saturation distribution between the VL and CL simulations, we calculate the difference in plume thickness at each (x, y) location at different times; see Figure 16. Red blocks represent locations where the CL Pc simulation exhibits a thicker plume and blue blocks represent locations where the VL Zero simulation is characterized by a thicker column.

The differences between the two end-member simulations is similar when increasing the number of injectors from one to four. The regions in the vicinity of the injectors display a much thicker plume in the CL Pc simulation, with a maximum difference of 160 m. At early times (Fig. 16a), there is faster lateral spreading plume in the VL Zero simulation, following the faster initial vertical migration. The cells further away from the injector are characterized by a large negative difference in plume thickness, which corresponds to the blue colour. As more CO₂ is injected, the regions nearer the wells generally have a thicker CL Pc plume, whereas the central part displays a thicker VL Zero plume. The absolute difference in plume thickness decreases at later stages (Fig. 16e, f), which reflects the progressive filling of the reservoir.

Overall, the results demonstrate that there is a noticeable difference in plume thickness between the VL Zero and CL Pc simulations. Therefore, disregarding the upscaled impact of capillary heterogeneity will likely result in erroneous saturation distribution predictions, with a significant under-prediction in plume thickness near the injectors.

This study has demonstrated the impact that capillary heterogeneity has on CO₂ plume migration in the dome structure and well configuration of the Endurance CO₂ storage site. We have adopted the well rates, well locations and reservoir properties from the publicly available Endurance reports as far as is possible.

The key qualitative impacts of the heterogeneities are to increase the swept plume volume and area. This is a result of intra-reservoir baffling and lateral preferential flow paths slowing the vertical segregation of the CO₂ from the in situ brine. This emphasizes the importance of including capillary heterogeneity upscaling in field scale simulations. Neglecting capillary heterogeneity results in very different plume dynamics. The impacts of these heterogeneities are pronounced for the earlier stages of injection, and decrease as the structure begins to fill. There are also differences in gravity current migration rates at the caprock but these are largely mitigated by the dome shape of the site, and would likely fall within the uncertainty range of forecast and history match simulations using conventional approaches.

These results and those discussed in the SI also indicate that the leading edges of the plume and the thin preferential flow pathways that form within the reservoir will likely fall under the seismic detection limit. In particular a significant amount of the plume (reaching 1400 m length for the single-injector scenario) may be undetected. These features are important for estimating the plume extent and migration speed, ultimately to determine any containment risks. Therefore, it will be important to utilize a range of monitoring techniques to track the plume in time.

Whilst a change in injector location results in variations in the total plume footprint area and swept volume, the overarching impact of capillary heterogeneity remains. We observe a similar impact of capillary heterogeneity after increasing the number of injectors from one to four. The five-well field development plan was designed partially to provide ‘better mitigation against any unforeseen field heterogeneities’ (BEIS 2021b). However, our results show that the impact of capillary heterogeneity does not diminish in this well configuration. Instead, a larger number of preferential flow pathways form and the individual plumes interact, which increases the complexity of the plume shape.

Our results show the significant differences in plume shape and speed arise from disregarding the upscaled impact of capillary heterogeneity on the multiphase flow functions in the geological setting and well configuration of the Endurance CO₂ storage site. Our work clearly demonstrates that it is crucial to properly account for capillary heterogeneity and its impact on multiphase flow in a full 3D reservoir model in order to understand and accurately model CO₂ migration.

Multi-layer plume formation has been observed at other CCS sites around the world, similar to that seen in our simulations for Endurance. For example, the plume at the Sleipner site in Norway has been described as having a ‘Christmas Tree’ shape (Chadwick and Noy 2010) consisting of nine separate layers characterized by high CO₂ saturations. Such multi-layer plumes may result in enhanced capillary trapping and dissolution due to the increase in plume area. This has important implications for CCS operations as injection could be designed to maximize preferential flow pathway formation and thereby improve storage efficiency.

The effect of upscaling capillary heterogeneity in Endurance is to preferentially increase the speed of the CO₂–water front along the bedding plane compared to the cross-bedding. This is because the pseudo relative permeabilities for the parallel to bedding plane direction do not result in a shock front from CO₂ displacing water, whilst the cross-bedding plane pseudos do result in a shock. This also means that the along bedding plane CO₂–water displacement is more unstable and prone to viscous fingering and the formation of preferential flow paths.

Disregarding the upscaled impact of capillary heterogeneity will likely result in under-prediction of migration speeds and spreading. Unexpected fast migration speeds have been observed at CCS sites worldwide. In Salah is one example, where the plume was found to have three times faster migration rates than had been predicted previously with flow simulations (Ringrose et al. 2009).

Based on the Endurance simulations, the leading edges of the plume will likely fall below the seismic detection limit (see SI) due to the thin, laterally extensive preferential flow pathways that form as a result of capillary heterogeneity. This is an important finding given that this leading edge of the plume also moves more rapidly than is predicted without consideration of capillary heterogeneity. Seismic surveys will likely produce erroneous CO₂ location estimations, in particular they will not necessarily identify the extent of the plume and thus fail to provide an early warning of it encroaching on leaky wells, faults or spill points. However, we also observe that capillary heterogeneity reduces the total number of locations which exhibit plume thicknesses smaller than the seismic detection limit after 25 years of injection. Thus, the formation of a multi-layered plume, due to the presence of capillary heterogeneity, may in fact reduce the challenges associated with plume detection at later times.

These insights are summarized in a schematic in Figure 1. The maroon plume shape represents the conventional view of CO₂ migration away from an injector, where the plume gradually thins towards the gravity tongue (Krevor et al. 2015). Our results (pink plume) indicate that CO₂ plume migration is more complex under the influence of capillary heterogeneity. Near the injector, preferential flow paths are created, which improve sweep efficiency and increases capillary trapping. As we move away from the injector, the gravity tongue thins and we observe faster migration speeds.

The heterogeneities we accounted for at the Endurance site have a much larger correlation length parallel to the bedding planes (700–1000 m) than perpendicular to the bedding planes (4 m). In addition the bedding planes are parallel to the direction of flow and the reservoir dip is small. The effect of capillary pressure in this case is to reduce the tendency of the injected CO₂ to segregate from the water due to buoyancy. This is seen in the lower upscaled gas relative permeability in the vertical compared to the horizontal direction. Including smaller scale features of the sedimentology, e.g. mm-cm scale bedding textures, would likely amplify the affects we have observed in our simulations.

Of course, different reservoirs have been formed in different depositional environments with different styles of heterogeneity, and there will be a commensurate variation in the impact of capillary heterogeneity. In different depositional environments there will be different correlation lengths and degrees of anisotropy in correlation length for the metre scale heterogeneity. The bedding planes may have a different orientation to the principal flow direction. The impact of capillary heterogeneity will thus vary. In the end member of a fining upward sequence the capillary heterogeneity will further reduce segregation whilst in a fining downwards sequence segregation will be further enhanced. In fractured (dual porosity) systems it is likely that capillary heterogeneity will both increase segregation and result in faster plume migration. While capillary forces are absent from the fractures, the capillary contrast at the matrix-fracture boundary will make it even harder for the CO₂ to access the low permeability, waterwet matrix between the fractures.

In this work, we have extended a 2D capillary-limit, steady-state upscaling scheme to model the 3D upscaled impact of capillary heterogeneity on the multiphase flow properties, capillary pressure and relative permeability. Using the workflow, we have investigated the impact of capillary heterogeneity on plume migration at the proposed Endurance storage site.

The focus of our study was the impact of heterogeneity at an industrial scale geological carbon storage site using realistic injection rates and multi-point CO₂ injection. By means of an extensive core characterization study, we evaluated the properties of the injection interval at the Endurance site. This in turn allowed us to model the likely flow behaviour in the reservoir. We investigated variations in flow behaviour with changes in well location and topography by replicating the proposed field development plan.

We find that the impact of capillary heterogeneity on plume migration differs between 2D and 3D models, primarily because the interactions of heterogeneities in all three dimensions affect the CO₂ flow pattern. This is particularly evident for the lateral migration speed, which is slowed down in 2D. These are novel insights that have previously been overlooked and clearly demonstrate the value of building a full 3D reservoir model.

Our work has shown that capillary heterogeneity has a significant impact on plume migration. Near the injector, it enhances the formation of CO₂ preferential flow pathways, which significantly improves vertical sweep efficiency and may result in an increase in capillary trapping and dissolution, if these effects were modelled. At distances far from the injectors, the gravity tongue is thinner and the plume travels faster as a result of capillary heterogeneity, although in this case the dome structure largely mitigated this speed enhancement. These are novel insights and have important implications for CCS operations. On the one hand, it may be possible to design well placement to take advantage of capillary heterogeneity to improve the swept volumes of the reservoir. On the other hand, it can lead to large variations in plume migration speed and needs to be properly accounted for in reservoir models.

We acknowledge Schlumberger and the Computer Modelling Group (CMG) for providing access to ECLIPSE and IMEX, respectively.

NW: conceptualization (equal), data curation (lead), formal analysis (lead), investigation (lead), methodology (lead), validation (lead), visualization (lead), writing – original draft (lead), writing – review & editing (lead); AM: conceptualization (equal), formal analysis (equal), methodology (equal), supervision (equal), validation (equal), writing – original draft (supporting), writing – review & editing (supporting); SJ: conceptualization (equal), formal analysis (equal), investigation (equal), methodology (equal), supervision (equal), writing – original draft (supporting), writing – review & editing (supporting); SA: formal analysis (supporting), investigation (supporting), methodology (supporting), writing – review & editing (supporting); SK: conceptualization (lead), funding acquisition (lead), resources (lead), supervision (lead), writing – original draft (supporting), writing – review & editing (supporting).

This work was funded through an iCASE studentship from the Engineering and Physical Sciences Research Council and bp (19000079), the UK Department for Energy Security and Net Zero (GB) (CCUS1014), and a Senior Research Fellowship from the Royal Academy of Engineering (RCSRF2223-1677).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The data associated with this work may be obtained from the following data repositories: https://doi.org/10.5285/7b616c86-e426-436d-87cd-1a3a8cc3c06d.