Modeling Pathways and Stages of CO 2 Storage

The storage of supercritical CO2 in deep geological formations can be partitioned in three stages: diffusion, early and late convection. The mass transfer of gaseous CO2 into an underlying brine layer is determined by the relative relevance of interacting multiphysics processes at the interface and the build up of a boundary layer. Within the brine layer convection emerges as a multiphysics phenomenon of coupled flow and transport in porous media. For the characterization of the three stages we use numerical experiments with perturbations of a reference homogeneous situation. We explore the effect of different type and size of perturbations. The simulations show that the onset of the convection state depends strongly not only on the perturbations, but also on settings of the numerical method. Moreover it is found that the early convection state may consist of several peaks and is thus more complex than in the idealized simple concept of a single peak. For the late convection stage the decrease of the total mass transfer into the system is generally confirmed.


INTRODUCTION
Carbon dioxide sequestration in deep geological formations has been suggested as a way of reducing atmospheric CO2 concentrations [1].In a study examining the future energy scenarios of 16 major countries the Sustainable Development Solutions Network (SDSN) andthe Institute for Sustainable Development and International Relations (IDDRI) conclude that deep decarbonisation pathways (DDP) 'are needed for increasing the ambition of country commitments to reduce their greenhouse gas emissions under the United Nations Framework Convention on Climate Change (UNFCCC)' [2], i.e. to make the 'transition to a low-carbon economy consistent with the internationally agreed goal of limiting anthropogenic warming to less than 2 degrees Celsius (°C)'.
The concept is to inject supercritical CO2 into a deep sub-surface formation that is filled with brine.The mass transfer of gaseous CO2 across the interface into the brine includes gas dissolution and diffusive transport.The brine in contact with CO2 will soon be saturated with CO2 and due to the diffusion of gas, a boundary layer will be formed that grows with time [3].The density of the water-CO2 solution increases with increasing CO2 concentration.
A second stage follows when the diffusion boundary layer becomes unstable.Convective motions appear as a multiphysics phenomenon of coupled flow and transport in porous media.Dense brine is moving downward and complementary lighter fluid is moving upward [4].The thickness of the appearing diffusion layer in the different stages and the concentration profiles within are affected by the flow pattern in the brine formation.In that ___________________________________ *Corresponding Author: ekkehard.holzbecher@gutech.edu.omway natural convection enhances the mass transfer of CO2 into the formation substantially.In this stage natural convection is the dominant process governing the mass transfer of CO2 into the brine.
Experimental and numerical studies have shown that the convective mixing appears mainly in large fingers [5].Downward moving fingers and transients with increased density reach the base of the permeable formation.In the third stage the convection eddies become weaker, as the density differences in the entire system decrease.With increasing carbon trapped in the subsurface brine reservoir, the convective motions weaken, causing the mass transfer across the interface to decrease.Finally, diffusion will again become the dominant mechanism.
Hassanzadeh et al. [6] characterize the three major stages during the process of carbon storage as follows: (A) diffusive mixing, (B) early convective mixing, (C) late convective mixing.
Figure 1.Stages in carbon storage: (A) diffusion, (B) early convection, (C) late convection; modified from: Hassanzadeh et al [6] Some characteristics of the three stages are sketched in Figure 1.Most crucial is the dimensionless Sherwood number Sh. Sh represents the total CO2 mass transfer into the brine layer.In the Sh number the mass transfer is normalized to diffusive mass transfer into a fully developed diffusion layer.A rigorous definition of Sh is given below.
In stage (A) the mass transfer due to diffusion is significantly higher than 1, as the concentration gradient is high in the early development of a diffusion boundary layer.We may distinguish an initial part of stage (A), which is only governed by diffusion and a second part, in which diffusion is still dominant, while convective movements in the system are in a state of early appearance.In the phase of early convection (B) Sh increases significantly.After reaching a peak, in the third stage (C) the mass transfer decreases gradually due to weakening convection motions.
The general behaviour of the entire dissolution process can be simulated by numerical modelling.In our model approach we use COMSOL Multiphysics to solve a system of partial differential equations, which describes coupled flow and transport in porous media.A derivation of this formulation is given in chapter 2. In chapter 3 we describe the model set-up and numerical details.In chapter 4 results are presented with focus on the different stages of the entire dissolution and storage of CO2 in the subsurface formation.

DIFFERENTIAL EQUATIONS
In variable density problems flow and transport are linked by a two-way coupling -in contrast to usual flow and transport problems with only a one-way link [7].The general flow equation is given by: with porosity , density , viscosity , pressure p, the gravity vector g and the permeability tensor k.In porous media the relevant Darcy-velocity v is given by the equation: The mass transport equation is described by: where c denotes CO2 concentration and D the general dispersion tensor.Transport is always dependent on the solution of the flow equation, as the velocity v determines the advective part.The link from transport back to the flow equation, which can be neglected in usual applications, is given through the dependencies (c), the dependency of fluid density on CO2 concentration.
The entire system of equations ( 1) to (3) is simplified as follows.Instead of the full dispersion tensor we use an effective diffusion coefficient of the fluid/solid system D. Another common simplification is the so called Oberbeck-Boussinesq assumption, which states that density differences are relevant at first place in the so called buoyancy term of the Darcy equation and can be neglected everywhere else.Thus the transport equation can be simplified to: For general analyses it is convenient to make a transformation into dimensionless form.
In 2D the flow equation can be formulated in terms of the streamfunction  (instead of pressure): where x and z denote the horizontal and vertical space directions, k is the corresponding permeability [7,14].Velocities and streamfunction are connected by the formulae: With the transformations , and correspondingly one obtains instead of equation ( 4) and instead of equation ( 5): A linear state equation for normalized concentration is assumed: Assuming constant permeabilities and viscosities equation (8) transforms to: The final step is the normalization of space ,which leads to: The dimensionless coefficient on the right hand side is the Rayleigh number.It is named after Lord Rayleigh, who published a paper on convection in free fluids in 1916, in which he introduced a dimensionless combination of parameters as characterizing of the convection problem in free fluids [8].The analogous dimensionless combination for the porous medium Ra, as denoted in eq. ( 11), was derived by Lapwood in 1948 [9].
Altogether the transformations from the original to the final dimensionless formulation, are: The given or very similar descriptions are discussed in the literature under the terms 'density-driven flow' [10,11,14], 'variable density flow and transport' [12] or 'natural convection in porous media' [11,13].In case of high Rayleigh numbers a very similar description is derived for 'density fingering' [14].Similar multiphysics formulations have been used dealing with other density related problems in aquifers, as saltwater intrusion [15,16], salination from saltdomes [17], thermal convection, etc [14,18,19,20].
Here we focus on combined mass transport and fluid fluxes, only.However and other aspects may also be taken into account in addition.Geomechanical processes are dealt with by Bornara et al. [21], Aker et al. [22].Thermal effects of heat transport are neglected in the derivation above.For the consideration of temperature effects in the borehole see Sponagle et al. [23].Double diffusive effects, arising from coupled mass and heat transport, are studied by Bouzgarrou et al. [24] and Islam et al. [25].The effects of anisotropic dispersion are examined by Ghesmat et al. [26], as well as by Musuuza et al. [27,28].In this paper inhomogenities are taken into account, as they trigger convective motions.A more general approach of heterogeneous layers is presented by Farajzadeh et al. [29].

Conceptual Model
In our model the problem of mixing is simulated in a 2D vertical cross-section through the porous medium formation below the interface between the overlying CO2 layer.For simplicity it is assumed that the interface at the top of the model region is horizontal, which also holds for the base of the modelled porous formation that is underlain by an impermeable geological barrier.Moreover, the porous medium itself in the model region is assumed to be homogeneous.These settings are common in studies on fingering effects of CO2 sequestration.
Above the interface a CO2 mass fraction of 0.0493 is assumed, following Pau et al. [30].The following parameters are taken from the same reference, adopted from the Carrizo-Wilcox aquifer in Texas.
Combining the physical parameters as shown in equation ( 11) for an assumed aquifer thickness of 100 m, results to Ra=4307.The Rayleigh number may change by several orders of magnitude, mainly due to permeability that may differ by orders of magnitude from one geological layer to another.The dependence of CO2 sequestration on the Rayleigh number is extensively discussed by Hassanzadeh et al. [31].Here we demonstrate and examine other features of the numerical simulations for an example Rayleigh number of Ra=5000, a value that has also been selected for numerical simulations by other authors.Connected with the parameters of Table 1 this corresponds to an aquifer thickness of 116 m.As model region we choose a representative cut out of permeable layer of the deep subsurface.As the emerging flow patterns are slender, i.e. are small in horizontal direction relative to the height, we use a square as model region, which captures the fingering phenomena that are representative for the entire formation.A sketch of the cross-section is depicted in Figure 2. The streamfunction approach has the advantage that the boundary conditions become simple: one can use the Dirichlet condition at all edges, if the model region is closed at all sides.There is no fluid flux across the boundaries; they are streamlines.In the model we use the flow boundary condition   0 everywhere.
For transport we chose a Dirichlet condition for CO2 concentration at top.This is more realistic than a flux boundary condition, as the fluxes of CO2 across the interface are not known beforehand.Mass fluxes depend on the concentration distribution in the diffusion layer below the interface, which is also influenced by the convective motions.Thus fluxes can be expected to change spatially and temporally with the flow pattern, and are an output of the model approach.
Except from the top boundary there are 'no flow' conditions for transport at all other edges, i.e. .At the vertical sides the mentioned boundary conditions for flow and transport reflect periodic conditions.At the bottom boundary they represent the interface to an impermeable layer, across which there are no fluxes, neither of fluid, nor of CO2.
For the onset of convection small disturbances of homogeneous conditions are of immense influence.As a nucleus these disturbances trigger convection.For that reason several types of disturbances are explored here.As convective motions set in near the top interface, disturbances in the diffusion boundary layer can be expected to be most relevant.We examine disturbances in boundary conditions at the interface boundary (section 4.1) and in initial conditions in the diffusion layer (section 4.2).Permanent disturbances in the entire domain are mentioned briefly (section 4.3).
The trivial initial conditions in the model are: zero CO2 concentration and no flow in the entire model region.Oscillatory and white noise disturbances of the initial concentration distribution are examined, which trigger convective motions.Another option is to use a boundary layer solution as initial condition for c, which is also studied.
Disturbed initial conditions reflect changes of local conditions at start of the simulation.Inhomogeneities of parameters, which may cause such disturbances of initial c are not explicitly taken into account in that case, and their effect on the further development can thus not be studied in that way.Thus it is surely more realistic, to introduce disturbances in parameters that remain present during the simulation.In our models disturbing the Ra number, appearing in equation, does this (11).

Numerical Model
For the numerical simulations we use the COMSOL Multiphysics [32] code, a versatile and flexible code for Finite Element solutions of coupled partial differential equations.The software can be applied to all kinds of (multi)-physical settings.It is equipped with a graphical user interface that allows easy handling and coupling of different physics modes.
For the mathematical formulation the dimensionless formulation of equations ( 7) and ( 11) is taken.We choose the Poisson equation as a classical differential equation for the flow equation (11), and the 'Convection-Diffusion' mode for the transport equation (7).The coupling variables from flow to transport are the velocity components.The transport to flow coupling is given by the buoyancy term on the right hand side of equation (11).
For high Rayleigh numbers, as they emerge in the case of carbon sequestration, the flow pattern is unsteady, showing fluctuating and transient features.Their dependence on numerical features, as meshing, is explored here.
Concerning finite elements we use irregular triangular meshes and Lagrange elements of 2nd order.The dependence on mesh refinements is tested using three meshes of different mean mesh spacing.In the following the three different meshes are labeled as coarse, medium and fine.Some details of the meshes are listed in Table 2.The solution of the coupled system is obtained using the direct MUMPS solver with row preordering for the monolithic system that includes flow and transport in a single matrix (see Holzbecher [33]).
The maximum dimensionless simulation time is chosen as t=0.01.With the input parameters, mentioned earlier, this corresponds to a time of more than 210 4 years.The numerical simulation proceeds in time steps.These are calculated automatically in order to achieve a required accuracy of the solutions.At initial time and at onset of convection the increased dynamics of the system is reflected by small timesteps, while they increase towards the end of the simulation.

RESULTS
It is well known and demonstrated in our simulations that light disturbances near the CO2brine interface grow into prominent fingers.These fingers subsequently merge to form larger fingers, developing into complex convective flow patterns, as shown in the sequel.While this holds for the real physical phenomena, it is also true for the numerical models.We discuss several ways in which disturbances can be introduced in the simulation runs and examine their influence on the resulting flow patterns.
We examine disturbances in boundary conditions, in initial conditions and in the domain.In the models all types of disturbances can be introduced simultaneously.However, we use a single type of disturbance only in each model run.Concerning the disturbed initial state, we consider oscillatory and random disturbances.

Boundary Disturbances
As boundary condition at the upper boundary, i.e. at the interface between the brine formation and the injected CO2, in the models we require: The first summand represents the elevated concentration; the second describes disturbances of oscillatory nature that are assumed to remain constant during the simulation.ε is the amplitude of the disturbance, which is assumed to be small.N represents the frequency of the spatial oscillations.As reference the values ε = 1/1000, N = 10 are chosen.
The disturbance is consistent with the side boundary conditions / 0 c n    .In some variations we also used the sin-function (instead of the cos in formula ( 13)), which is not consistent with the side boundary conditions.The inconsistent disturbances cause stronger fingers to appear at the vertical edges (not shown here).
Disturbances at the boundary may reflect inhomogeneities of the porous medium (porosity, permeability, diffusivity) at the interface.They persist during simulation time.However, the parameter variations do not enter the model directly.Equation ( 13) reflects the effect of these inhomogeneities on concentration.How inhomogenities in the model domain can be considered in the model approach, is outlined in section 4.3 below.
The different plots of Figure 3 show flow patterns shortly after start of the simulation.They are obtained for two different oscillatory disturbances at the boundary and for the three meshes, noted above (Table 2).
As Figure 3 shows the oscillatory boundary conditions is reflected by the number of fingers in the flow pattern, independent of the mesh size.However, deviances due to the mesh can be observed in details of the finger shape.The coarse mesh result for the longer wavelength (N=10) disturbance (top left) looks most irregular.For the shorter wavelength (N=18) the result (top right) looks more regular, but every second period is amplified differently.For the shorter wavelength the flow patterns of medium and fine mesh (center and bottom right) are almost identical.But for the longer wavelength there are still recognizable differences in the CO2 concentration distributions within the fingers: the medium mesh result (center left) shows fingers, in which the high concentration region is split into an upper and a lower part.The fingers seem to be near to release a transient heavy bubble to sink.Such a split is not visible in the fine mesh result (bottom left).

Initial Disturbances
It is well known that initial disturbances play an important role concerning the onset of convection, not only in the physical reality but also in numerical modeling.Using COMSOL Multiphysics Holzbecher [34] highlights this for the advance of a moving front of a fluid with different viscosity.Here similar effects can be observed also for a fluid with changing density.We study initial disturbances of two types: oscillatory, similar to the boundary disturbances, and random white noise.Moreover we study two types of boundary layer representations.

Oscillatory Disturbances 4.2.1.1. Steep Gradient
As in the previous section oscillatory disturbances are specified by an amplitude (=maximum size) and a frequency, ε and N: The second factor in formula ( 14) describes an initial steep concentration gradient in vertical direction below the interface.Not very far away from the interface the concentration and disturbance become almost zero.A different physically motivated representation of the boundary layer will be discussed in the following sub-section.
Figure 4 shows the development of convection starting from initial disturbances, at two time instants (t=0.0015 and 0.0023) and the three meshes described above.Colours represent the distribution of CO2.Arrows indicate flow size and direction at some locations.The c=0.1 contour line is drawn in white colour.The figure shows how the small fingers near the CO2-brine interface (see Figure 3) grow into prominent fingers and transients, sinking towards the base of the aquifer.
In the upper two figures (computed with the coarse mesh) an irregular pattern of fingers can be observed.From the earlier state (left) to the later state (right) a merging of fingers seems to appear.Three adjacent fingers of similar size arrange themselves to a pattern, showing only three major fingers.Such behaviour was also observed by Gesmat et al. [26].In this simulation the c=0.1 contour has not reached the base yet.
The results for the same model run with the medium mesh is shown in the center figures.The finger pattern is obviously much more irregular.From earlier to the later state transients form, i.e. plumes of increased concentration, which then loose the connection with the finger, from which they emerge.At later time (center right), the first transient already hits the base with a concentration higher than 10% CO2.A second transient is likely to reach the bottom soon.
In the shown simulations the c=0.1 contour arrives at the bottom at dimensionless time t=0.0024 for the coarse mesh and at t=0.0021 for the medium mesh.For the boundary condition disturbances, as discussed above, and initial white noise disturbances (see below) we found only slightly deviating values.These times match with the arrival time calculated by Farajzadeh et al. [35] for the same Rayleigh number.
The output for the simulation on the finest mesh for the same physical set-up and disturbances, is depicted in the lower two sub-plots.The finger pattern looks very different from the coarse and medium mesh calculations.Here the oscillatory features are mainly kept.Identical fingers develop and move towards the base.Between the two time instants shown the fingers become more pronounced in the upper part.In the lower part the c=0.1 contour is far from reaching the base; in fact only a slight movement of its mean position can be seen in that time period.
Figure 5 shows the further development of the reference case for two further time instants (t=0.0035 and 0.005), again showing also the mesh dependencies.More CO2 has been dissolved leading to increased CO2 concentrations.The distributions for the coarse and the medium mesh show only few spots with normalized concentrations lower than 0.1 (top & center left) for the earlier time.At t=0.005 the minimum concentration is 0.3 in the coarse mesh (top right), and 0.23 in the medium mesh (center right) simulation.Regular flow features can not be observed in these simulations.
The different flowpath for the simulation with the fine mesh that was already seen in the previous figure is traced further in Figure 5 (bottom).The concentrations are much below those of the other two simulations, due to the obviously weaker convection of the finger pattern.For the earlier time the c=0.1 contour has not reached the aquifer base yet.At t=0.035 a first deviation from the regular finger pattern can be observed, near to the left and right boundaries.These irregularities increase.At t=0.05 the first finger has already hit the base, and a second is to follow soon after.At this time instance stronger convective patterns emerge.
Altogether two pathways can be observed.The simulations with the coarse and medium mesh show similar results with a quick breakdown of the initial flow pattern and transition to a strong convection regime with high mass transfer.In the fine mesh simulation a regular finger pattern is preserved for a much longer time.However, finally there also is a transition to strong irregular convection features.
where y denotes the dimensionless space variable.The concentration distribution according to equation ( 15) can be taken as initial condition for the simulation.In our model we combine this with disturbances, i.e we start the simulation with a disturbed diffusion layer.
For oscillatory disturbances the initial condition thus becomes: where t0 denotes the assumed length of the a first period with pure diffusive regime, i.e. the initial part of stage (A) (see Figure 1).As in the early time of simulation diffusion remains the only relevant process, i.e. the entire length of phase (A) is longer than t0.The length of this second part of stage (A), that lasts until the transition from diffusive to convective regime, i.e the beginning of stage (B) is dependent on the parameter t0.This is explored below.Figure 6 depicts the development of CO2 transfer across the top interface.The top figure results are obtained with the coarse mesh, the bottom figure with the medium mesh.Time t is counted from the end of the purely diffusive stage.The mass transfer is represented by the Sherwood number Sh, given by integrating the normal fluxes, in the dimensionless system given by: The transition from diffusive to convective regime is obvious in all graphs of the figure.During the initial stage (A) Sh is slightly decreasing and rises sharply with the onset of convective movements.While within stage (A) we see Sh below 20, a rise to values above 100 can be observed in the early convection stage (B).The order of magnitude for the simulated duration of the diffusion phase is in the order or 10 -4 dimensionless time units.Obviously the initial diffusive mass transfer and also the length of diffusion dominance, i.e. of stage (A) depend on the parameter t0: the lower t0 , the higher the peak mass transfer.
The length of stage (A) can be obtained by adding the duration of the purely diffusive regime t0 and the simulated time until onset of convection.With increasing t0 the simulated duration until onset of convection is also increasing.Physically this is a reasonable behavior: with longer t0 there is a higher mass transfer into the system, which causes the concentrations gradients to decrease with time, and smaller gradients are unfavorable for the onset of convective motions.For the prediction of the entire length of stage (A) this an inconvenient result, as the length of the first part of the pure diffusive stage t0 is not known.
Figure 6 also shows a mesh dependency.The medium mesh results have two time periods with rising Sh, in contrast to the coarse mesh results.Peak Sh also depends on the mesh.Not mesh dependant are the Sh levels in the stage of declining convective motions (C), and these are also independent of t0.In the figure this is clear for the coarse mesh results, but holds also for medium and fine mesh simulations for longer simulation times.Moreover, for fixed parameter t0 the onset of early convection stage is not mesh dependent, which is also confirmed by runs with the fine mesh (not shown).
Figure 6: Fluxes across the interface, measured by the Sherwood number Sh, in dependence of t for varying boundary layer build-up time t0; top for coarse mesh, bottom for medium mesh, calculated for oscillatory initial condition 4.2.2.Random Disturbances Real geological systems are never homogeneous, and disturbances of homogeneity are most likely not oscillatory.In the numerical model we may introduce usually unknown disturbances by random distributions.Following Riaz et al. [36] we utilize initial white noise in the boundary layer for this purpose.In the diffusion layer, the cos term in the equation ( 16) is replaced by a local random number.
For the three simulations of the reference set-up with different meshes Figure 7 depicts the mass transfer changing with time.The fluxes are represented by the Sherwood number Sh, as introduced by equation ( 17).For comparison the curves for the three runs of the previous sub-section are added, which were started with oscillatory disturbances according to formula (14).Three stages, diffusion (A), early convection (B) and late convection (C) can be identified for most of the simulation runs.However, the stages do not follow the simple shape sketched in Figure 1.Simulated stage (A) is relatively short.Its length obviously depends on the prescribed initial condition.For the runs with oscillatory disturbances it lasts longer than for random disturbances in the boundary layer.The reason for the deviances is not the type of disturbance, but the approach concerning the diffusion layer.As outlined, in the shown results oscillatory disturbances are obtained with the steep gradient formulation, whereas random disturbances result from the diffusion layer initial condition.
For the runs with random disturbances the length of the diffusion stage is the same.The runs with oscillatory disturbances show a complex mesh dependency: shortest duration of stage (A) is obtained for the coarse mesh, the longest for the medium mesh.
Concerning the convection stages the simulation results deviate from the ideal simple concept of a branch with increasing fluxes, followed by a branch of decreasing fluxes, as shown in Figure 1.In fact the fluxes fluctuate significantly.That may be explained due to the dynamic behaviour of the fingering pattern, in which single fingers play an important rule for a certain time period before they are diffused.These fingering phenomena influence the flow pattern and thus also the mass transfer across the boundary.
There is a strong increase of the boundary flux with the onset of convection: in all simulations Sherwood numbers rise above 100.The decreasing trend in the final period of the simulations can also be observed without doubt.Thus the general existence of stages, as proposed and shown in Figure 1, can be confirmed.But the transition between early and late convection is not easy to localize.After the first maximum the simulations often show another increase of mean Sh.The fluxes stay at a high level and in most simulations drop at dimensionless time t=0.003.The range of Sh in the late convective state is neither influenced by the assumptions concerning the diffusion layer, nor by the mesh size.
One simulation clearly follows a different pattern.It is the fine mesh simulation that starts with oscillatory disturbances in the diffusion layer.Obviously this represents a different pathway in stages (A) and (B).The flow patterns of that pathway are already visualized in Figures 4 and 5: the bottom figures there show that the regular finger pattern, introduced by the oscillatory boundary conditions, remains stable for a much longer time.

Domain Disturbances
The random disturbances, reflecting inhomogeneities in geological characteristics, were introduced for the initial condition so far.In that way they form the nucleus of convective flow patterns, but have no further influence after the start.In real geological systems instabilities may arise foremost from inhomogenities in geological characteristics.Geological inhomogenities, as for example deviations from the mean hydraulic permeability, remain constantly present during the entire time-period.In the numerical models such inhomogeneities can be represented by disturbances of flow or transport parameters.Within the analytical approach, presented in chapter 2, this is best mimicked by introducing a perturbation to the Rayleigh-number Ra.In the model such disturbances are then permanent.
We ran several simulations with white noise disturbances of Ra with disturbances size given by a parameter ε similar to equations ( 14) and ( 16).For small value of ε=1/1000 as default, as in formerly presented perturbations, the qualitative results were not different from the ones reported so far.The examination of effects from larger size perturbations, as they are even likely in geological systems, can be subject of further studies.

CONCLUSION
High Rayleigh-number coupled flow and transport phenomena, as they can be expected for CO2 storage in deep sub-surface geological strata, show a complex dynamics.Our studies focused on the stages of pure diffusion (A), early convection (B) and late convection (C), as sketched in Figure 1.The mass transfer depends on the flow patterns with fingers and transients that emerge from differing local conditions.In mathematical models such differences are represented by disturbances of homogeneous conditions.By numerical simulations we examined several types of disturbances in order to study their influence on the different stages of the storage process.Disturbances can be introduced in boundary conditions, initial conditions or they can be permanent within the domain.
Duration of length (A) is usually in the order of magnitude of 10 -4 dimensionless time units.Duration of stage (B) is one order of magnitude longer than stage (A).There is a transition to stage (C) at approximately 310 -3 dimensionless units.Stage (C) is clearly longer than stage (B).At t=10 -2 the CO2 concentration has risen above 60% of saturation even at the base of the geological formation.
Assumptions concerning the initial development of the diffusion layer have an effect on the critical time appearing in the simulation, marking the transition from diffusive to convective regime.Strong gradients as they appear initially favor the onset of convective instabilities.If there is longer time for the diffusion layer to develop, the gradients are smaller, which also leads to a longer critical time within the simulation for the start of the convective regime to become dominant.From this we conclude that the length of stage (A) can not be predicted by the model, if the diffusion layer equation (15) with unknown t = t0 is applied.In contrast, model runs with the steep gradient approach deliver a length of stage (A) independent from mesh size.
The early convective stage also depends on the assumptions concerning the build-up of the diffusion layer.In the early convective mixing the Sherwood numbers increase up to above 100 (for Ra=5000).Some runs show that mean values of mass transfer often remain at a plateau until decrease, others show a second higher peak.The simulation results deviate from the ideal simple concept of a branch with increasing fluxes, followed by a branch of decreasing fluxes, as shown in Figure 1.However, the general existence of three stages can be confirmed, but the situation is more complex than in the ideal simple sketch.
Two pathways could be identified.They differ mainly concerning the stability of the initially developed regular finger pattern.For the run with finest mesh and oscillatory initial conditions we found that the regular pattern with fingers of almost equal size (length and width) remains stable for a much longer time.In all other runs single fingers, that develop earlier, arrive earlier at the base and induce the entire system to develop strong convective circulation patterns.For the two pathways the stages differ substantially!Although the differing pathway is obtained with the finest mesh, i.e takes smaller scale details into account, we assume that the resulting pathway is most likely unrealistic.Geological media are generally characterized by strong inhomogenities.They are much bigger than the small random fluctuations used here.The pathway with early onset of convection seems to by more realistic in that respect.Physically the pathway with single fingers is more likely, as hydrogeological inhomogenities will surely prevent the emergence of regularly fluctuating fronts.Out of the presented runs only one showed an exceptional pathway.Most other simulation runs, which are not reported here in detail, also show the early instability of the initial pattern.Only few runs follow the exceptional pathway, whose strictly regular pattern is surely artificial.

Figure 2 .
Figure 2. Schematic view of coordinate system, model region and boundary conditions

Figure 7 .
Figure 7. Flux across interface (Sherwood number) in dependence of normalized time; results for reference case for different meshes and initial conditions (see text)

Table 1 .
Characteristic values of relevant parameters for CO2 sequestration

Table 2 .
Characteristics of finite element meshes and resulting degrees of freedom for the coupled model, using quadratic elements