Some thoughts on Darcy-type ﬂ ow simulation for modelling underground CO 2 storage, based on the Sleipner CO 2 storage operation

We take three ﬂ ow simulators, all based on Darcy ’ s Law but with di ﬀ erent numerical solver implementations, to assess some of the issues surrounding their use to model underground CO 2 storage. We focus on the Sleipner CO 2 injection project, which, with its seismic monitoring datasets, provides unique insights into CO 2 plume devel- opment during a large-scale injection operation. The case studies ﬁ rstly compare simulator performance in terms of outputs and run-times on carefully matched model scenarios; then we compare numerical with analytical Darcy solutions to explore the potential for modelling simpli ﬁ cation; ﬁ nally we look at the e ﬀ ects of including conservation of energy in the simulations. The initial case-study used simpli ﬁ ed axisymmetric model geometry to simulate the upward ﬂ ux of CO 2 through a heterogeneous reservoir, incorporating multiphase ﬂ ow with coupled CO 2 dissolution into formation brine. All three codes produced near-identical results with respect to CO 2 mi- gration velocity and total upward CO 2 ﬂ ux at the reservoir top. The second case-study involved 3D modelling of the growth of the topmost layer of CO 2 trapped and migrating beneath topseal topography. Again the three codes showed excellent agreement. In the third case-study the simulators were tested against a simpli ﬁ ed analytical solution for gravity currents to model the spreading of a single CO 2 layer beneath a ﬂ at caprock. Neglecting capillary e ﬀ ects, the numerical models showed similar layer migration and geometry to the analytical model, but it was necessary to minimise the e ﬀ ects of numerical dispersion by adopting very ﬁ ne cell thicknesses. The ﬁ nal case-study was designed to test the non-isothermal e ﬀ ects of injecting CO 2 into a reservoir at non-ambient temperature. Only two of the simulators solve for conservation of energy, but both showed a near identical thermal anomaly, dominated by Joule-Thomson e ﬀ ects. These can be signi ﬁ cant, particularly where reservoir conditions are close to the critical point for CO 2 where property variations can signi ﬁ cantly a ﬀ ect plume mo- bility and also seismic response. In conclusion, the three simulators show robust consistency, any di ﬀ erences far less than would result from geological parameter uncertainty and limitations of model resolution. In this respect the three implementations are signi ﬁ cantly di ﬀ erent in terms of computing resource requirement and it is clear that approaches with simpli ﬁ ed physics will pay rich dividends in allowing more detailed reservoir hetero- geneity to be included. Contrary to this, including conservation of energy is heavier on computing time but is likely to be required for storage scenarios where the injectant stream is signi ﬁ cantly di ﬀ erent in temperature to the reservoir and most critically for shallower storage reservoirs where CO 2 is close to its critical point.


Introduction
Numerical simulations are a vital tool for understanding the short, medium and long-term fate of CO 2 injected in the subsurface. Indeed, one of the key regulatory requirements outlined in the European Directive on CO 2 storage (EU, 2009) is to demonstrate "…. conformity of the actual behaviour of the injected CO 2 with the modelled behaviour". It is important therefore to determine the validity and applicability of numerical simulators for reproducing the key coupled processes related to flow of CO 2 in the reservoir.
Here we look at three commonly used numerical simulators to assess their efficacy for modelling underground CO 2 storage. The assessment is in three parts. First we test the comparative performance of the simulators on two identical model scenariosone with axisymmetric geometry, one with full 3D geometry. Second we compare the numerical results with analytical solutions for a spreading CO 2 layer. Finally we assess the degree to which it is necessary to include thermal effects in CO 2 flow simulation.
Various authors (e.g. Pruess, 2005;Class et al., 2009) have run code comparisons on CO 2 injection problems; some have focussed on hypothetical test cases and others such as Singh et al. (2010) have utilised calibration data from a real injection project. In all these cases a satisfactory match between the codes was obtained but with some differences that can largely be ascribed to variation in modelling implementation between the operators at the various institutes carrying out the code comparison.
Our comparison differs somewhat from previous studies in that all three simulators were tested in-house (at BGS), with model meshes carefully designed to be as identical and/or equivalent as practicable. In this way, modelling differences due to 'operator-effect' were largely eliminated.
For this study we deployed three simulators: TOUGH2 from Lawrence Berkeley National Laboratory (Pruess, 2004); ECLIPSE 100 the industry-standard black oil simulator from Schlumberger (Schlumberger 2011); PFLOTRAN an open-source parallel subsurface flow and reactive transport code (Lichtner et al., 2015). All three are based on established multi-phase implementations of Darcy's Law for fluid flow in a porous medium.
The CO 2 storage project at Sleipner (Baklid et al., 1996) has a comprehensive time-lapse seismic monitoring programme which provides unique detail on large-scale CO 2 plume development. The casestudy simulations presented here are based around these monitoring results, for overall realism, and to enable us to explore some of the issues around the numerical simulation of actual physical processes. But it is important to stress that we do not attempt to obtain exact historymatches of the Sleipner monitoring data; this is a multi-faceted task outside the scope and aims of our comparative study.
1.1. The Sleipner CO 2 storage operation CO 2 separated from natural gas produced at the Sleipner field in the North Sea (Norwegian block 15/9) is being injected into the Utsira Sand, a regional saline aquifer of late Cenozoic age (Fig. 1). The aquifer comprises mostly unconsolidated sand of high porosity (> 30%) and high permeability (> 1 Darcy) and is generally in excess of 200 m thick in the Sleipner area. A number of thin intra-reservoir mudstones, typically 1-2 m thick, are evident from geophysical logs acquired in wells around Sleipner (Fig. 1).
The CO 2 is injected via a deviated well in a dense phase at a depth of 1012 m below sea level, approximately 200 m below the top of the aquifer. Injection commenced in 1996 at a roughly constant rate, with around 15 million tons of CO 2 stored by 2016. The growth of the resulting plume has been monitored using repeat time-lapse 3D seismic data acquired in 1999, 2001, 2002, 2004, 2006, 2008 and 2010 (Fig. 2).
The CO 2 is imaged on the seismic data as a number of high amplitude sub-horizontal reflections within the reservoir (Fig. 2). Most of this reflectivity is thought to arise from thin layers of CO 2 trapped beneath the intra-reservoir mudstones which are partially but not wholly sealing (Arts et al., 2008;Chadwick et al., 2004Chadwick et al., , 2005Chadwick et al., , 2010Chadwick and Noy, 2015;Chadwick et al., 2016). The seismic data suggest that CO 2 has migrated vertically upwards through the reservoir via a chimney in the mudstones (Fig. 2) located a little to the south of the injection point. Nine interpreted reflective layers had formed by 1999 (when CO 2 first reached the top of the reservoir), and each individual reflective layer has remained mappable on all of the subsequent surveys.

Simulators used in the code comparison
ECLIPSE 100 is a fully implicit, 3D, 3-phase, isothermal black oil simulator. Fluid Pressure-Volume-Temperature (PVT) properties are interpolated from look-up tables as a function of pressure. The simulations in this code comparison were run using block-centred geometry (Schlumberger, 2011), with inter-block transmissibility calculated using an arithmetic average of the cell interface area coupled with the harmonic average of the permeability. The density and compressibility of CO 2 as a function of temperature were calculated using the Span and Wagner (1996) equation-of-state and converted to the relevant black oil representation using the scheme published by Hassanzadeh et al. (2008): brine is represented as an oil phase in the simulator, allowing for CO 2 dissolution. Activity coefficients for CO 2 and H 2 O were obtained directly from the equation-of-state using the methodology described by Spycher and Pruess (2005), while the solubility of CO 2 in brine was calculated according to Duan et al. (2006). Brine density and viscosity were taken from the International Association for the Properties of Water and Steam tables, using Ezrokhi's method to calculate the density effect of salt and dissolved CO 2 (Zaytsev and Aseyev, 1992). The viscosity of CO 2 was calculated as a function of temperature and density using relationships published by Vesovic et al. (1990) and Fenghour et al. (1998).
TOUGH2 was designed as a general purpose, multi-phase, non-isothermal simulator (Pruess et al., 1999;Pruess, 2004). It uses integral finite differences to achieve spatial discretisation, with fully implicit time stepping. By default it implements an upstream weighting scheme to calculate inter-node transmissibility coefficients, but this was changed to a harmonic averaging scheme to match the ECLIPSE models. The code comparison used the ECO2N fluid properties module (Pruess, 2005), which uses the methodology described by Spycher and Pruess (2005) to compute activity coefficients. CO 2 density and viscosity are derived from correlations published by Altunin (1975); which provide a very close approximation to the Span and Wagner (1996) equation-ofstate. These PVT correlations are used to generate tables of CO 2 and brine properties that are then interpolated during the computation to obtain specific phase parameters. The solubility of CO 2 in brine is calculated according to Duan et al. (2006).
PFLOTRAN is an open-source, parallel subsurface flow and reactive transport code, built on top of the PETSc family of PDE numerical showing reservoir heterogeneity (γ-ray logs on the left tracks and resistivity logs on the right tracks). The reservoir sand has characteristically low γ-ray and resistivity readings so peaks within the sand denote thin mudstones.
solvers (Lichtner et al., 2015). It contains a specific multiphase model of the brine − CO 2 fluid system (MPHASE), with CO 2 density calculated using the Span and Wagner (1996) equation-of-state. MPHASE calculates the viscosity of CO 2 as a function of temperature and density using relationships published by Fenghour et al. (1998). Again, the solubility of CO 2 in brine is calculated according to Duan et al. (2006). PFLO-TRAN uses the International Formulation Committee of the Sixth International Conference on Properties of Steam (Haywood, 1965) to compute the density and viscosity of water. These PVT correlations are used to generate tables of CO 2 and brine properties that are then interpolated during the computation to obtain specific phase parameters.

Reservoir and fluid properties
Reservoir temperatures were based on robust long-term measurements from the nearby Volve water production operation, quoted in Alnes et al. (2011). Reservoir pressures were assumed to be initially hydrostatic. Internal fluid property modules supplied with TOUGH2 and PFLOTRAN were used to compute the density and viscosity of CO 2 as a function of temperature and pressure (see above), with example values shown in Table 1. Brine fluid properties were calculated for a sodium chloride mass fraction of 0.032, giving a mean brine density and viscosity of approximately 1020 kg m −3 and 8 × 10 −4 Pa.s respectively. The black oil parameters used in ECLIPSE100 were calculated using CO 2 densities from the Span and Wagner (1996) equationof-state (see Table 1) using the scheme described by Hassanzadeh et al. (2008).

The case-study simulations
Four case-study simulations have been used for our code comparison, each of which was based around aspects of the Sleipner CO 2 injection operation described above. The first case-study aimed to approximate whole plume geometry in a vertically heterogeneous reservoir using a 2D radial axisymmetric model. The second case-study investigated the temporal evolution of the topmost CO 2 layer in threedimensions (see Chadwick and Noy (2010) for more information). This model is similar to the 'Sleipner Benchmark' released by Statoil in 2011 and described by Singh et al. (2010) and Cavanagh (2013). In the third case-study the TOUGH2 simulator was compared with an analytical solution for a single CO 2 layer spreading beneath a flat caprock. The final case-study, a simple axisymmetric model of a warm buoyant CO 2 plume rising upward through a uniform sandstone reservoir, was used to test the thermal modelling capabilities of two of the codes (PFLO-TRAN and TOUGH2).

Model setup
The axisymmetric model for the full Sleipner plume used here has previously been described by Chadwick et al. (2010) and Chadwick and Noy (2015). It was created for the TOUGH2 simulator and incorporates a 2D radial axisymmetric mesh, with 80 cells in the X (I) grid direction increasing in width logarithmically with distance from the injection well. The mesh is only one cell thick in the Y (J) direction, but the cell volume increases with X so that it represents the full circumference at the radial distance X from the injection well.
The simulation grid is divided vertically into 170 cells with variable dimensions, chosen to give detailed resolution in the vicinity of the injection point and to include the thin mudstone layers known to trap CO 2 in the reservoir. The top of the reservoir is positioned at a depth of 797 m below sea level, with the base of the reservoir at 1083 m. Noflow boundary conditions were placed at the top, base and perimeter of the model, although the grid extends to a radial distance of 2.5 × 10 5 m, effectively behaving as an "open aquifer" over the injection timescale.
The geological model comprises 16 stratigraphical layers: four 'lower' sands and five 'upper' sands, with intervening thin mudstones and the topmost sand unit capped by 50 m of very low permeability caprock (Table 2). Sand permeabilities are consistent with measurements of 1.6-3.3 Darcy from a single core  but Fig. 2. Cross-section from the time-lapse 3D seismic data at Sleipner showing the baseline (1994) dataset and selected repeat surveys. The amplitudes have been normalised using the maximum survey amplitude to facilitate comparison between vintages. Note strong reflections corresponding to the CO 2 plume (top panels). Maps of the plume expressed as normalised average absolute reflection amplitude in a travel-time window containing the plume (bottom panels). Line of cross-section is marked with a black line in map view. Location of the plume feeder chimney is indicated by a black circle in map view and a black arrow on the seismic sections. Black polygons from 2001 onwards mark the extent of the topmost CO 2 layer within the overall plume footprint.

Table 1
Representative CO 2 fluid properties used in the code comparison. The values were calculated using the Span and Wagner (1996) equation-of-state. Note that TOUGH2 uses correlations published by Altunin (1975) to obtain CO 2 fluid properties, however the two formulations show excellent agreement within the PT range of the case-studies. These properties were converted into black oil tables for ECLIPSE using the scheme described by Hassanzadeh et al. (2008).

Depth
Temperature ( modified following Chadwick and Noy (2015). The permeability of the intra-reservoir mudstones was adjusted to allow CO 2 to reach the top of the reservoir by 1999, matching time-lapse seismic observations which showed that CO 2 had reached the caprock just prior to the (October) 1999 repeat survey. This constraint has resulted in a much higher effective permeability for the intra-reservoir mudstones than has been measured on intact core samples of Sleipner caprock (around 4.0 × 10 −7 Darcy, Harrington et al., 2010;Chadwick and Noy, 2015). This 'semi-permeable' behaviour might reflect preferential flow through pervasive small-scale fractures resulting in very low effective threshold pressures (Cavanagh and Haszeldine, 2014) or higher permeability feeder chimneys in the mudstones (Hermanrud et al., 2009). Relative permeability and capillary pressure curves for CO 2 and brine in the Utsira Sand were computed by fitting a Van-Genuchten model to measurements from core samples measured in the laboratory (Eric Lindeberg personal communication). A similar function was used for the mudstones, but with an order-of-magnitude increase in capillary entry pressure (see Table 2 and Fig. 3).
A variable injection rate with a mean value of 27 kg s −1 was used, based on the actual values measured at the wellhead (Ola Eiken, Statoil, personal communication, see Fig. 3d).

Results
Sample results from the axisymmetric model comparison are presented as CO 2 saturation profiles extracted at simulation time-steps corresponding to the 1999, 2001 and 2006 seismic monitor surveys (Fig. 4). For comparison with the observed data, markers show the equivalent radial extent of each seismically mapped CO 2 layer, calculated by fitting a circle with the same area as the observed layer area. It is clear that the different codes produce remarkably consistent CO 2 distributions in the plume, with individual layer spreading very consistent between the models. The models also provide a reasonable match to the layer extents mapped on seismic data but it is noticeable that with time the modelled radii tend to progressively outstrip the observed equivalent radii. This is to be expected, given that layer spreading in the axisymmetric model is beneath perfectly flat mudstone seals, whereas the real seals have a degree of undulation, with buoyant ponding tending to retard spatial spread.
The relative consistency of each simulator in terms of advective fluid mass transfer through the plume can be tested using the volumetric growth of the topmost CO 2 layer and it is clear that the codes all show excellent agreement in predicting this (Fig. 5). CO 2 volume in the topmost layer can also be derived from time-lapse seismic analysis (see Chadwick and Noy, 2015 for methodology) and it is clear that the model results are consistent with seismic observations, with the proviso that very thin CO 2 accumulations at low saturation (< 0.2) cannot be reliably mapped on the time-lapse seismics (to account for this, only CO 2 saturations > 0.2 have been included in the simulated layer volume curves).

Model setup
Recent quantitative seismic analysis of CO 2 plume migration at Sleipner has focused on the topmost spreading layer in the plume (e.g. Chadwick and Noy, 2015;Williams and Chadwick, 2012;Furre et al., 2015). This layer lies directly beneath the reservoir caprock and is clearly and stably imaged on time-lapse 3D seismic data, which provides reliable constraints for flow simulation.
The monitoring data show that growth of the topmost layer between 1999 and 2006 was accomplished by rapid lateral spreading of CO 2 to infill a topography of domes and ridges beneath the reservoir topseal (Fig. 6).
The geological model for the code comparison is based on that developed by Chadwick and Noy (2015) and comprises a single uniform sand unit, 16 m thick, with a porosity of 0.37 and a permeability of 8 Darcy. Relative permeability and capillary pressure curves for CO 2 and brine in the Utsira Sand were computed by fitting a Van-Genuchten model to measurements from core samples measured in the laboratory (Table 2 and Fig. 3a). The assigned permeability is substantially higher than core-plug measurements (∼3 Darcy), but recent re-assessment of regional permeability and wireline log data suggest that the permeability of the topmost sand body might be significantly higher than previously supposed (Williams and Chadwick, 2017). This is examined further in the Discussion.
The top reservoir surface was mapped from the baseline (1994) seismic dataset and depth-converted using a layer-cake velocity model derived from available well data. The model was discretised using 60 50 × 50 m cells along the X axis and 111 50 × 50 m cells along the Y axis. The 16 m thick reservoir was divided vertically into 8 × 2 m cells. No-flow boundary conditions were placed at the top and base of this layer (simulating an impermeable caprock and underlying mudstone layer), while the lateral model domains were maintained at near hydrostatic pressure conditions. This was achieved by using large boundary cells (TOUGH2 and PFLOTRAN) or a pore-volume multiplier (ECLIPSE).
The flow model assumed a single feeder to the top layer of CO 2 positioned at the prominent gas chimney observed on seismic data (Fig. 2), with CO 2 flux matched to volumes derived from the seismic data (Fig. 5). Temperature in the sand unit was set to 29C in the simulators and fluid properties were again calculated using the Span andWagner (1996) equation-of-state in ECLIPSE andPFLOTRAN andAltunin's (1975) correlations in TOUGH2 (Table 1).

Results
Results from the three models show the extent of the top spreading layer at the time of the 2006 time-lapse seismic survey (Fig. 7). All models show the dual effects of radial spreading from the CO 2 source/ feeder with evidence of buoyancy-driven migration in the north-east where the CO 2 is starting to migrate northwards beneath the linear ridge (Fig. 6).
As before, all three codes show good agreement in CO 2 distributions but with some variation in local migration patterns at the layer edges. In the Eclipse 100 simulation more CO 2 has spilled beyond the CO 2water contact at the eastern margin of the grid compared to the TOUGH2 or PFLOTRAN models. This reflects the fact that Eclipse 100 is a black oil simulator that assumes a constant reservoir temperature. This was set to the average temperature in the layer (∼30C), whereas a temperature gradient was used in the TOUGH2 and PFLOTRAN model runs. All of the simulation models are consistent with each other, but in absolute terms they fail to match the observed extents of layer spreading; this aspect is examined further in the Discussion.

Case-study 3: comparison of numerical codes with an analytical solution
The third case-study compares the numerical simulations with an analytical solution employing similar but simplified physics for the spreading of a single CO 2 layer beneath a flat caprock. Lyle et al. (2005) and Bickle et al. (2007) presented analytical solutions for gravity flows in a permeable medium with axisymmetric symmetry. Their model (Fig. 8) comprises a porous permeable medium filled with a fluid into which a less dense fluid is introduced along a line source under an impermeable caprock. The fluid ponds under the caprock and spreads radially. To obtain the solution several simplifying assumptions were made, including no capillary pressure, no viscosity difference, no relative permeability effects and constant fluid densities.

Analytical solution
Darcy's Law and the continuity (conservation of mass) equation give key relationships between the rate of fluid input (αQt (α−1) , where a value α = 1 corresponds to release at constant flux Q (m 3 /s)), the radius of the spreading layer, rN(t) (m), as a function of time t (s), and the thickness of the layer h (m) as a function of radial distance r (m), and time (s).
The radius of the layer is given by: And the thickness is given by: Where: Ø is the porosity, k is permeability (m 2 ), ρ is the density of the introduced fluid (kg/m 3 ), μ the viscosity of the introduced fluid (Pa s) and g' = gΔρ/ρ, the reduced gravity with Δρ the difference in density between the introduced fluid and the initial fluid filling the medium (kg/ m 3 ). The similarity variable η is given by: η N (α) (i.e. η at r = r N ) is a function of α only and is given by: Where the scaled similarity variable y = η/η N ranges from 0 at the source to 1 at the outer edge of the CO 2 layer and f(y) is given by numerical solution of the differential equation:  This is a good approximation for f except when y is small, and η N is well approximated by:  The equations predict, that 1) for constant input flux (α = 1) the radius of the CO 2 layer will be proportional to the square root of time; 2) for constant flux the axial thickness of the layer h(0,t) is nearly invariant; 3) there are simple relationships between both the radius and thickness of the layer and the input flux and physical parameters porosity, viscosity and permeability.

Numerical model setup
The numerical simulation of a single CO 2 layer was carried out in TOUGH2. An axisymmetric model was setup injecting dense-phase CO 2 into the base of an aquifer 112 m thick with an impermeable caprock. The mesh is highly refined, with elements 5 cm thick at the top of the aquifer and 1 m wide out to a radius of 250 m, then gradually expanding to the outer boundary at 20 km.
Key properties of the model are based broadly on Sleipner (Table 3), with fluid properties as follows: brine salt mass fraction 0.032; brine density 1020 kg m −3 ; brine viscosity 8.60 × 10 −4 Pa s; CO 2 density 720 kg m 3 ; CO 2 viscosity 5.93 × 10 −5 Pa s. The injection point was placed 102 m below the reservoir top with an injection rate of 10 kt per year, the latter being of the same order as the CO 2 supply to the topmost layer at Sleipner.
The analytical model assumes identical reservoir and fluid properties and conditions as the numerical model, with a constant flux of 10 kt per year of CO 2 along an axial line source (Fig. 9).

Results
In the TOUGH2 model the CO 2 rises in a thin buoyant column to the base of the caprock and then spreads radially beneath it (Fig. 9). The modelled layer is around 3-4 m thick in the axial part (above the feeder column) thinning gradually towards its leading-edge. A steepening of the CO 2 -water contact (CWC) towards the perimeter is an effect of capillary resistance to flow of the CO 2 , defined by the capillary pressure curve. CO 2 saturations in the layer, again defined by the capillary pressure curve, are generally high, but decrease towards the CWC across a zone termed the capillary fringe. The analytical model (shown as a black curve in Fig. 9) shows similar layer development and geometry to the numerical model but with significant differences. The CWC maintains a roughly uniform dip out to the layer edge and the layer spreads rather further than in the numerical model. Full CO 2 saturation throughout the layer is implicit in the analytical model, because it does not include capillary effects. Volumetrics require therefore that the layer in the analytical model be somewhat thinner than in the numerical model. The analytical solutions predict that for a constant injection rate the axial thickness of the layer is invariant (Eq. (2)). This behaviour is supported by the numerical model which maintains roughly constant axial layer thicknesses through time (Fig. 9).
To examine the effects of capillary forces on layer geometry and saturation a second scenario was developed, this time injecting 500 t of CO 2 per year (Fig. 10). A lower injection rate was selected because capillary effects are relatively more significant in thinner layers and the reduced injection rate provides this.
The first numerical model run included capillary pressure and shows a familiar geometry with the CWC steepening towards the layer edge with a saturation fringe (Fig. 10a). It is notable that the axial layer thickness is smaller than those above, reflecting the smaller input flux. The corresponding analytical model shows the familiar thinner layer with somewhat wider spread and uniform CWC dip (shown as a black line in Fig. 10a). The second numerical model run omitted capillary pressure effects from the computation (Fig. 10b). The resulting layer has full CO 2 saturation throughout and shows a much improved match to the analytical model in terms of both the layer spread and layer   Bickle et al., 2007).

Table 3
Physical properties used in the numericalanalytical comparison.
In order to test for the effects of numerical dispersion, the model mesh was modified such that the cell heights at the top of the reservoir were increased from 5 cm to 25 cm (Fig. 11). Compared to the detailed numerical model this had the effect of increasing the volume of each modified grid cell and increasing the minimum layer thickness. One consequence is to reduce the free phase CO 2 saturation in the layer (Fig. 11b) because the same amount of CO 2 is distributed within a greater volume in each cell. Another consequence is to reduce the layer lateral spread at the top of the plume. This is an important point because accurate modelling of lateral plume migration is a key requirement for robust long-term storage prediction. It is notable that even in the modified model (Fig. 11b) the layers are still very thin compared with cell dimensions typically adopted for generalised reservoir flow modelling, so this numerical dispersion effect is likely to be significant in many flow models.
To further the numerical/analytical comparison, repeat model runs were carried out in Eclipse 100 and PFLOTRAN with results that were, to all intents and purposes, identical to those from TOUGH2. It is evident therefore that, at least for the growth of a single layer, the numerical and analytical approaches give very comparable results. In the numerical models, capillary pressure effects clearly influenced layer geometries and saturations, especially for thinner layers, but when these were removed (to match the analytical model assumptions), results were almost identical.

Case-study 4: thermal effects
Much of the flow simulation work traditionally carried out for CCS assumes, for simplicity, that CO 2 is injected at ambient reservoir conditions. This is not usually the case. Depending on wellhead temperature and specific adiabatic conditions in the wellbore, CO 2 will generally not have the same temperature as the reservoir at the injection perforations. So, for example, when injecting into strongly depleted fields, adiabatic decompression in the wellbore will lead to injecting CO 2 significantly colder than the reservoir. In contrast, at Sleipner adiabatic compression in the wellbore leads to the injection of CO 2 warmer than the surrounding reservoir. The temperature of the CO 2 at the injection perforations is estimated at ∼48C or slightly above (Alnes et al., 2011), around 13C warmer than ambient reservoir temperature at the injection point. Furthermore, interpretation of the time-lapse gravity data (Alnes et al., 2011) suggest that the average density of the CO 2 plume at Sleipner is compatible with its having a warm, less dense, axial core.
The effects of a plume of warm buoyant CO 2 rising to the top of the reservoir without cooling to ambient reservoir temperature are worthy of investigation. Two of the codes compared in this study (TOUGH2 and PFLOTRAN) solve for conservation of energy and are suitable for modelling heat propagation in a migrating CO 2 plume. A new axisymmetric model was set up, injecting CO 2 at a rate of 27 kg s −1 , but at 49C, into the base of a homogeneous sandstone reservoir 275 m thick and overlain by 750 m of low permeability shale caprock (Table 4). This is simpler than the model shown in Fig. 4 in that there are no intra-reservoir mudstones to interrupt the plume ascent to the reservoir top.
The CO 2 saturation distributions produced by the two simulators in this simple reservoir model are to all intents and purposes identical (Fig. 12), a simple buoyant axial column of CO 2 rising above the injection point and feeding a single CO 2 layer spreading radially beneath the caprock.
Fluid temperature distributions from the two models are also very Fig. 10. Growth of a single CO 2 layer with an injection rate of 500 ty −1 a) Numerical model with capillary pressure b) Numerical model without capillary pressure. The black line in a, b shows the CO 2 -water contact predicted by the analytical solution. Fig. 11. The effects of numerical dispersion on growth of a single CO 2 layer with an injection rate of 500 ty −1 a) CO 2 saturation after 6 years for a model with a mesh spacing of 100 × 5 cm at the top of the reservoir b) CO 2 saturation after 6 years for a model in which the grid cell height has been increased from 5 to 25 cm. similar, comprising an axial column of elevated fluid temperatures above the injection point and a radially spreading thermal anomaly at the reservoir top. It is noticeable that the thermal imprints of the rising column and the spreading layer are of lesser spatial extent than the corresponding fluid saturation features. This is because of radial heat loss from the column into the surrounding reservoir and vertical heat loss from the spreading layer into the caprock and underlying reservoir.
The thermal anomaly at the injection point is +14C, corresponding to the temperature difference between the injected CO 2 and the ambient reservoir. Above the injection point the temperature anomaly reduces to +8C at the reservoir top and reduces further radially, such that it becomes negligible at distances more than about one third of the full layer spread More detail of the temperature evolution is revealed by a vertical temperature profile through the grid cells immediately above the injection point at a distance of 2.5 m from the model axis (Fig. 13a). Both codes show a consistent temporal and spatial evolution of the thermal anomaly. The dominant physical processes occurring in the model are rapid vertical advection of warm CO 2 combined with cooling due to expansion as the hydrostatic pressure drops with elevation above the injection point. Steady-state heat transfer occurs after ∼365 days, and thereafter the vertical temperature profile closely approximates the Joule-Thomson (JT) cooling curve (Fig. 13a) with a coefficient of 5C/ MPa for the warm buoyant CO 2 (Span and Wagner, 1996).
The models show that cooling due to expansion is significant and takes place a lot faster than conductive heat transfer to the surrounding reservoir. The CO 2 has cooled to a temperature of ∼37C by the time it has reached the top of the reservoir (some 8C above the background reservoir temperature) and the thermal anomaly propagates laterally beneath the caprock for a distance of 700 m from the axis (Fig. 13b).
Because this particular model scenario is based on Sleipner which is a shallow reservoir with the injected CO 2 close to its Critical Point, the effect of the thermal anomaly on CO 2 properties (Fig. 14) is significant. A ∼50% decrease in viscosity and density in the core of the plume both contribute to increased mobility of the CO 2 phase. The bulk modulus (Fig. 14d) is also markedly reduced which will modify the seismic properties of the CO 2 . This will affect the rock physics of CO 2 saturated reservoir rock and should be taken into account for quantitative analysis of the seismic response such has been attempted a number of times at Sleipner (e.g. Chadwick and Arts, 2005;Furre et al., 2015).

Validity of Darcy-flow assumption
The three simulators show results that are remarkably consistent  with each other. Differences in the modelled outer edges of the layer extents are generally around 50 m or less, though locally up to ∼100 m. This is similar to or less than the uncertainty of layer edge detection on the seismic data itself, Bickle et al. (2007) noting that the outermost parts of the seismically-imaged layers would not be detectable due to their being too thin, with under-estimation of layer extents by up to 100 m. In absolute terms however none of the simulation casestudies match the observed monitoring data very well, notably in their failure to reproduce the rapid lateral spread of the layer observed by 2006. History-matching the observed growth of this layer has proved challenging generally, with most published simulations having difficulty in predicting the rapid northward migration of CO 2 along the prominent linear ridge in the base of the topseal (Figs. 6 and 7). Various authors have investigated the effects of small uncertainties in feeder geometry, topseal topography, permeability anisotropy and gas composition (e.g. Chadwick and Noy, 2015;Zhu et al., 2015) but the Darcybased modelling approaches have had mixed results in replicating the observed layer growth rates. Our numerical and analytical models are all based around a 'Darcy' type flow formulation. It has been proposed that because CO 2 migration in the topmost layer is dominated by buoyancy, rather than pressuredriven viscous forces, Darcy-type flow models, strictly valid only for slow viscous fluid flow, might not be wholly appropriate and modelling with alternative physics should be utilised (Cavanagh, 2013;Cavanagh and Haszeldine, 2014). In this context flow models with alternative physics have been proposed, notably the invasion-percolation (IP) scheme (Cavanagh, 2013). At Sleipner, IP modelling allows good spatial/geometric matching of the topmost layer spread, but it is a quasistatic technique that does not incorporate a temporal element (the fluid overcomes capillary resistance essentially instantaneously) and so is not well suited to the type of temporal history-matching that time-lapse monitoring requires (Oldenburg et al., 2015). Recent laboratory experiments (Krishnamurthy et al., 2017) comparing core-flood CO 2 saturation measurements against Darcy and IP simulations have found a tendency for IP to markedly underestimate CO 2 saturations, leading Krishnamurthy et al. (2017) to suggest hybrid modelling approaches, such as using multiple IP runs to pre-condition reservoir grid properties (e.g. capillary threshold pressure variation) prior to full Darcy simulation. More fundamental modelling approaches, reaching down into pore-scale processes, are also being developed, such as Lattice-Boltzmann modelling (e.g. Sukop and Thorne, 2007) but these are computationally very demanding and not currently suitable for simulating two-phase flows in reservoirs with complex geology and high resolution geological characterisation datasets.
Recent support for the Darcy-based modelling approach is provided by detailed re-assessments of the reservoir properties at Sleipner (Williams and Chadwick, 2017, Cowton pers. comm.) which indicate significant depositional heterogeneity in the top reservoir sand, with wide permeability variation. When these more heterogeneous reservoir properties are incorporated into Darcy flow models a much closer history-match can be obtained (Williams and Chadwick, 2017).
It is clear that computing power is a primary limit on simulation accuracy, notable in limiting the resolution and geological detail that can be incorporated into the reservoir meshes. Our comparison shows marked differences in processor requirements between different Darcy simulators, but it is also clear that Darcy modelling is susceptible to useful simplifications and the vertical equilibrium approximation for numerical modelling has proved useful in simulating the Sleipner plume (e.g. Nilsen et al., 2011). More radical is the development of purely analytical solutions to model CO 2 flow (e.g. Nordbotten et al., 2005;Lyle et al., 2005;Bickle et al., 2007). These allow extremely rapid computation of layer flows, albeit in a simple reservoir, and our comparison shows that published analytical solutions are consistent with a full physics Darcy numerical model. More recently a hybrid approach has been developed, whereby analytical solutions for gravity flows beneath dipping seals have been discretised within a finite-difference numerical scheme to allow very short computational times even on a PC (Laurence Cowton pers. comm.). This level of computational efficiency allows the incorporation of complex geology at the small (metres) length scales that are really required for realistic modelling.
There is therefore a crucial trade-off between including full complex physics in flow models and the need to include complex and high resolution geological heterogeneity. It seems that Darcy modelling with suitable simplifications to allow radically reduced computation times is the best way to proceed at the present time, particularly where geological uncertainty is significant and where the spread of a migrating CO 2 layer is accurately tracked through time.
The inclusion of thermal effects is clearly important, and is rarely included in CO 2 plume simulation. For injection into depleted gas fields CO 2 entering the reservoir will generally be significantly cooler than the surroundings and accurate modelling of the evolving CO 2 front will be essential. Even in the normally-pressured aquifer of the Utsira Sand, temperature has an important effect on fluid mobility Williams and Chadwick (2017), and as we show here, potentially on the seismic response of the fluid-rock system.  14. (a) Thermal anomaly, (b) CO 2 viscosity, (c) CO 2 density and (d) CO 2 bulk modulus resulting from CO 2 injection into a homogeneous sandstone reservoir. The output time is equivalent to the 2006 seismic monitor survey. The grid cell containing the injection point is marked with an arrow in the lower left corner of each plot. G.A. Williams et al. International Journal of Greenhouse Gas Control 68 (2018) 164-175 6.2. Model runtimes ECLIPSE 100 is at least 10 times faster than either PFLOTRAN or TOUGH2. This large difference in performance reflects the fact that both PFLOTRAN and TOUGH2 use numerical differentiation for the solution of the non-linear problem, while ECLIPSE uses analytical derivatives. This can bring significant computational saving. In addition TOUGH2 and PFOTRAN do not implement a linear solver optimised for the type of multi-phase problems encountered in reservoir simulations. ECLIPSE uses a nested factorisation that can offer significant speedup compared to standard linear solvers.

Conclusions
Three numerical multiphase flow simulators, one black-oil (ECLIPSE100) and two fully compositional (PFLOTRAN and TOUGH2), were tested on their ability to model CO 2 injection into a saline aquifer. Modelling case-studies were based on the Sleipner injection operation, which provided an appropriate context for large-scale storage, in terms of CO 2 layer dimensions, thickness profiles and spreading rates, plus excellent time-lapse monitoring data. All case-studies were set-up carefully on the three simulators to virtually eliminate 'operator effects'.
The simulators all showed excellent agreement in modelling the upward flux of CO 2 through an isothermal heterogeneous reservoir and in each case modelled saturation profiles showed good agreement with the distribution of CO 2 observed on time-lapse seismic surveys. 3D models of the topmost CO 2 layer, migrating beneath the caprock, and clearly imaged by seismic data, also showed excellent consistency between the three codes. Minor differences related to solver implementations, meshing, smoothing and small equation-of-state discrepancies. In all cases modelling differences were less than monitoring uncertainty. In absolute terms however none of the case-studies produced close matches of the layer's temporal evolution, notably its rapid lateral migration and very high mobility. One possibility is that the physics of fluid flow in this layer is different to the Darcy physics used by the simulators, but our favoured explanation is that imperfect reservoir characterisation is the root cause.
Numerical simulations of a simple, thin spreading layer with an analytical model also showed an excellent match and highlighted the effects of capillary forces in terms of layer fluid saturation and geometry. The importance of fine-scale modelling was also evident, as coarser numerical meshes resulted in strong numerical dispersion with impaired layer spreading accuracy.
Two of the codes (TOUGH2 and PFLOTRAN) were also run in nonisothermal mode to investigate the effects of injecting CO 2 significantly warmer than the ambient reservoir temperature. Again, the codes showed excellent agreement in predicting the distribution and magnitude of the resulting thermal anomaly and its temporal development. It is clear that, particularly in shallow reservoirs where stored CO 2 is close to its Critical Point, thermal effects are significant, because they affect CO 2 mobility and also its seismic properties, important for quantitative monitoring.