Solute Driven Transient Convection in Layered Porous Media

Abstract

CO₂ geological sequestration has been proposed as a climate change mitigation strategy that can contribute towards meeting the Paris Agreement. A key process on which successful injection of CO₂ into deep saline aquifer relies on is the dissolution of CO₂ in brine. CO₂ dissolution improves storage security and reduces risk of leakage by (i) removing the CO₂ from a highly mobile fluid phase and (ii) triggering gravity-induced convective instability which accelerates the downward migration of dissolved CO₂. Our understanding of CO₂ density-driven convection in geologic media is limited. Studies on transient convective instability are mostly in homogeneous systems or in systems with heterogeneity in the form of random permeability distribution or dispersed impermeable barriers. However, layering which exist naturally in sedimentary geological formations has not received much research attention on transient convection. Therefore, we investigate the role of layering on the onset time of convective instability and on the flow pattern beyond the onset time during CO₂ storage. We find that while layering has no significant effect on the onset time, it has an impact on the CO₂ flux. Our findings suggest that detailed reservoir characterisation is required to forecast the ability of a formation to sequester CO₂.

1 Introduction

The behaviour of CO₂ in a brine-saturated porous medium during CO₂ geological storage involves several processes. According to IPCC [1], a combination of traps safely store injected CO₂ within an underground storage formation. The vertical migration of buoyant CO₂ is restricted to lateral spreading by the formation seal (structural trap), the imbibing resident fluid can cut-off and immobilize trailing CO₂ plume within formation pores (residual trap), the dissolution of CO₂ in resident fluid removes CO₂ from a highly mobile phase (solubility trap), and the reaction of CO₂ saturated brine with formation rock securely stores CO₂ (mineralization trap). Tracking the buoyant migration, spreading, dissolution, and reaction of CO₂ simultaneously is non-trivial since they occur at different temporal and spatial scales. However, understanding the subsurface behaviour of CO₂ is crucial to controlling the global climate change through safe implementation of CO₂ geological sequestration.

To overcome the modelling challenges, many researchers focused on one or more of the trapping mechanisms at a time [2,3,4,5]. Our interest is on solubility trapping which progresses to trigger density driven convection. CO₂ dissolution rate, which measures solubility trapping, has been estimated for Sleipner field in Norway [6] while downwelling convective fingers are visually observed in experimental studies on convection in Hele-Shaw cells [7]. Previous studies of convection in homogeneous systems use the onset time when the instability commences and the critical wavelength of the convective fingers to characterise convection [2, 3]. The flow pattern after the onset time can be described in several regimes: diffusive, flux growth when the convective instabilities grow, constant flux when the fingers progress towards the bottom of the domain, and final period (when flux gradually declines) [4]. In heterogeneous systems, many studies on transient convection are in systems with random permeability field or impermeable dispersed barriers [5]. Layering which exist naturally in sedimentary rocks has not received worthy attention, but it can affect the dynamics of fluid flow due to the change in permeability. Therefore, we use numerical simulations with perturbations from numerical artefacts to show that layered systems can alter CO₂ flux beyond the onset time.

2 Governing Equations and Model

$$ \frac{{\partial {\text{u}}}}{{\partial {\text{x }}}} + \frac{{\partial {\text{v}}}}{{\partial {\text{z }}}} = 0 $$

(1.1)

$$ {\text{u}} = - \frac{\text{k}}{{{{\upmu }}}}\frac{{\partial {\text{p}}}}{{\partial {\text{x }}}} $$

(1.2)

This problem of density driven convection in a brine saturated porous medium is modelled in single phase comprising the flow of brine and the transport of CO₂ dissolved in brine. Each model aquifer is 2-dimensional in the cartesian coordinate system (x-z plane). The aquifer has a thickness H, impermeable (u = v = 0) and constant CO₂ concentration (c = C_s) at the top boundary (z = 0), and impermeable and zero CO₂ concentration flux $\left( {\frac{\text{dc}}{\text{dz }} = 0} \right)$ at the bottom boundary (z = H). Initially (t = 0), the domain is assumed to contain no CO₂ (c = 0). This initial and boundary conditions are common in the subject of convection in the context of the geological sequestration of CO₂ [2,3,4,5]. The system of equations for this problem is

$$ {\text{v}} = - \frac{\text{k}}{{{{up\mu }}}}\left( {\frac{{\partial {\text{p}}}}{{\partial {\text{z }}}} - \left( { 1 + {{\beta c}}} \right){{\rho }}_{0 } {\text{g}}} \right) $$

(1.3)

$$ \emptyset \frac{{\partial {\text{c}}}}{{\partial {\text{t }}}} + {\text{u}}\frac{{\partial {\text{c}}}}{{\partial {\text{x }}}} + {\text{v}}\frac{{\partial {\text{c}}}}{{\partial {\text{z }}}} = \emptyset {\text{D}}\left( {\frac{{\partial^{2} {\text{c}}}}{{\partial {\text{x}}^{2} }} + \frac{{\partial^{2} {\text{c}}}}{{\partial {\text{z}}^{2} }}} \right) $$

(1.4)

$$ {\text{v}} = - \frac{\text{k}}{{{{\upmu }}}}\left( {\frac{{\partial {\text{p}}}}{{\partial {\text{z }}}} - \left( { 1 + {{\beta c}}} \right){{\rho }}_{0 } {\text{g}}} \right) $$

(1.5)

where $c,$ $k,$ $D,$ $\emptyset$, $g$,$\mu$ and $p$ denote CO₂ dissolved concentration, permeability, diffusivity of CO₂ in brine saturated porous medium, porosity, gravity, viscosity and pressure respectively. The Darcy velocity components are $u$ in the horizontal and $v$ in the vertical direction. Boussinesq approximation is made so that the dependence of density (${{\rho }}$) on the concentration ($c$) which is assumed to be linear ${{\rho }} = {{\rho }}_{0 } \left( { 1 + {{\beta c}}} \right)$ is applied only in the buoyancy term $\left( {{{\rho g}}} \right)$ in (3), where ${{\rho }}_{0 }$ is the density of pure brine and ${{\beta }}$ is the expansion coefficient. The length, time, and concentration scale are H, the diffusive timescale $\frac{{H^{2} }}{D}$, and $C_{s}$ respectively. The velocity scale is $\frac{\emptyset D}{H}$, and the dimensionless pressure is expressed as $p^{\prime} = \frac{{\left( {p - \rho_{0 } gz} \right) }}{\emptyset \mu D}$. The main dimensionless number that describes the relative importance of diffusion and convection is the Rayleigh number, $Ra = \frac{{k\rho_{0} \beta C_{s} gH}}{\emptyset \mu D}$. The system of Eqs. (1.1–1.4) is made dimensionless, linearized, and then analysed using linear stability analysis (LSA) based on the quasi-steady state approximation (QSSA) needed due to the transient base profile. The procedure and validity of QSSA have been discussed previously [3].

3 Numerical Simulations

The governing equations and boundary conditions are numerically solved using Eclipse 300 which is a multi-phase, compositional reservoir simulator. Each model is 2 m × 1.01 m in the horizontal (x) and vertical (z) directions respectively discretized into 200 grids in x $\left( {i = 1,2 \ldots 200} \right)$ and 101 grids in z $\left( {j = 1,2 \ldots 101} \right)$ directions, where $i$ and $j$ are the spatial discretizations. The initial and boundary conditions are as described in Sect. 1.2. We rely on the numerical artefacts to trigger the instability.

4 Results and Discussion

The relationship $t_{0}^{'} = a{\text{Ra}}^{ - 2}$ between the onset time ($t_{o}^{'}$) and the Rayleigh number in the homogeneous systems is obtained both in the LSA and in the numerical simulations (Fig. 1.1a). This relationship and the constant a ≈ 56 from the LSA in this work are similar to that obtained in the previous study [3]. However, in comparison to the theoretical result, the numerical result for the constant a ≈ 807 is large due to weak perturbations from numerical artefacts which can delay the onset time. This can also be due to the difference in the definition of the onset time in the LSA and simulation (Fig. 1.1). Beyond the onset time, LSA breaks down, and we depend on the numerical simulations to observe the flow regimes. We estimate the average flux (F) using the rate of change of the total dissolved concentration (CO₂ liquid mole fraction) in the domain given by $F = \frac{1}{2 }\frac{d}{dt}\left( {\mathop \smallint \limits_{0}^{2} \mathop \smallint \limits_{0}^{1} cdxdz } \right)$ which is mathematically equivalent to the average flux from the top boundary $\left( { - \frac{1}{2 }\mathop \smallint \limits_{0}^{2} \frac{dc}{dz}dx} \right)$. The diffusive, flux growth and steady flux regimes that we find (Fig. 1.1b) have been previously observed [4]. The scaling agreement between our numerical simulation and LSA results for the onset time, and the similarity of our simulation results with previous numerical studies for the flow regimes after the onset time validates our numerical models. Hence, we extend this analysis to include heterogeneous cases using our simulation models.

Fig. 1.1

a LSA and numerical simulation results for the relationship between the onset time and Rayleigh number in homogeneous systems for Ra = 322.97, 430.44, 751.87, 1070.63. b Average flux vs time plot for Ra = 0 (diffusive system) and Ra = 430.44 also showing the profile for CO₂ liquid mole fraction (c’) at the top. The variation in Ra is solely due to permeability difference. Both plots are in dimensionless unit. The onset time in the LSA is the time when the perturbations growth rate become positive while the onset time in the simulation is when the standard deviation of the horizontal concentration profile across the domain increases

The heterogeneous models each has a depth, H, and comprises two homogeneous media with one overlaying the other. The thickness of each of the homogeneous media is half of the depth of the full system $\left( {\frac{1}{2 }H} \right)$. We investigate the effect of the layered heterogeneity on the onset time of convective instability and on the subsequent flow pattern. The nomenclature, TDiff - B430, interpreted as a heterogeneous system in which the top strata is a diffusive system (Ra = 0) and the bottom layer has Ra = 430.44 is adopted.

The onset of convective instability is not affected by the layering (Fig. 1.2). The top layer determines the onset time, and the bottom layer has no role in this setting. Similar to the onset time, the time of the flux growth depends largely on the properties of the top layer (Fig. 1.3). The time of flux growth in the heterogeneous systems with Ra = 430.44 at the top are similar and the amplitude of the growth can be approximated by that in the homogeneous system with Ra = 430.44. However, after the flux growth, the flux profile depends on the Ra of the bottom layer when dissolved CO₂ reaches this layer. The system with the purely diffusive bottom layer has a sharp flux decline (Fig. 1.3) due to the restriction of flow from the permeable top layer by the diffusive bottom layer (see concentration profile in Fig. 1.2). When a high Ra layer lies above a low Ra layer (T1071-B430 and T751- B430), the width of the convective fingers reduces when the fingers reaches the bottom layer reducing the flux (Fig. 1.2). The flux growth is controlled by the permeability in the top region, and during late times, the flux can be approximated by the flux in the homogeneous system Ra = 430.44. These results suggest that we can fairly estimate the flux in some layered heterogeneous systems with the flux in homogeneous media at certain times.

Fig. 1.2

Numerical simulation results for the relationship between the onset time and Rayleigh number in heterogeneous systems for Ra = 322.97, 430.44, 751.87, 1070.63 when Ra = 430.44 is at the top and at the bottom. At the top of this figure, the profile of CO_2 mol fraction in liquid phase in T430-BDiff, TDiff-B430, T430-B1071 and T1071-B430 at ${\text{t}}^{\prime} = 0.018$ are presented

Fig. 1.3

Average flux vs time plot for the homogeneous case Ra = 430.44 and the presented heterogeneous cases

5 Conclusions

The role of layered permeability heterogeneity on CO₂ solute convection in a brine saturated geological porous medium has previously not received much attention. We developed a numerical model on Eclipse 300 and validated our numerical results against theoretical work for the onset of convective instability and against previous studies for the flow regime beyond the onset time. Our results indicate that the bottom layer has no significant effect on the onset of convective instability while dissolved CO₂ remains in the top layer. Beyond the onset time, the presence of the bottom layer affects the flow regimes. This result implies detailed knowledge of potential storage formations is required to successfully implement CO₂ storage.