Convective Mixing Driven by Non-monotonic Density

Abstract

CO$_\text {2}$ injection for enhanced oil recovery (CO$_\text {2}$-EOR) or for storage in depleted oil and gas reservoirs can be a means for disposing of anthropogenic CO$_\text {2}$ emissions to mitigate climate change. Fluid flow and mixing of CO$_\text {2}$ and hydrocarbons in such systems are governed by the underlying physics and thermodynamics. Gravity effects such as gravity override and convection are mechanisms that can alter fluid flow dynamics, impacting CO$_\text {2}$ migration, oil production and eventual CO$_\text {2}$ storage at the field scale. This study focuses on convection in a miscible setting caused by non-monotonicity in oil density when mixed with CO$_\text {2}$, i.e., a maximum mixture density occurs at an intermediate CO$_\text {2}$ concentration. We perform high-resolution simulations to quantify the convective behavior in a simple box system where gravity effects are isolated. We show that convection of CO$_\text {2}$ in oil is dependent on whether CO$_\text {2}$ originates from above or below the oil zone. From above, convection follows classic convective mixing but is accelerated by viscosity decrease with increasing CO$_\text {2}$. From below, convection flows upward due to CO$_\text {2}$ buoyancy, but is countered by downward convection due to the heavier mixture density. This convective system is significantly more complex and efficient than from above. We characterize the instabilities in both early- and late-time regimes and quantify mixing rates. For a 100 mD reservoir, convective fingers would be on the order of centimeters in width and mix over a meter length scale within days to a month, depending on the placement of CO$_\text {2}$. The simulations are performed in non-dimensional form and thus can be rescaled to a different reservoir parameters. Our results give important insights into field-scale impacts of convective mixing and can guide future work in development of upscaled models and experimental design.

1 Introduction

Storage of carbon dioxide (CO$_\text {2}$) in geological formations is a means to reduce atmospheric emissions of greenhouse gases. CO$_\text {2}$ capture and storage (CCS) with utilization for enhanced oil recovery, or CO$_\text {2}$-EOR, is perceived as a cost-effective method of disposing captured CO$_\text {2}$ emissions (Heidug et al. 2015). CO$_\text {2}$-EOR has been performed for many decades with a focus on hydrocarbon recovery (Hadlow 1992; Perera et al. 2016). However, with CCS, several factors will bring additional challenges to a well-established industry (Eide et al. 2019). First, CO$_\text {2}$ storage will be increasingly emphasized in order to meet greenhouse gas emissions targets, requiring operators to monitor and predict the long-term fate of CO$_\text {2}$ in the reservoir. Second, CO$_\text {2}$-EOR will be introduced into the offshore environment, where economic constraints may restrict the number and placement of wells thus impacting flow regimes in the reservoir. In both cases, the interaction between CO$_\text {2}$ and hydrocarbons at the fine scale plays an important role, and therefore, detailed understanding is required for effectively managing CO$_\text {2}$ storage efficiency in CO$_\text {2}$-EOR reservoirs, or in any storage reservoir where existing hydrocarbons are present but not produced.

CO$_\text {2}$-EOR is most commonly applied as a miscible flood, but can also be applied in an immiscible setting. In miscible flooding, the injected CO$_\text {2}$ mixes with in situ hydrocarbons into a single phase (Holm 1986). The absence of interfacial tension allows CO$_\text {2}$ to mobilize residual oil, creating an oil bank that is driven toward producers (Lake 1989). Viscosity reduction and oil swelling are the main drivers in a miscible CO$_\text {2}$ flood (Green and Willhite 1998). Density changes also occur within the mixing zone, but is considered a secondary effect in viscosity-dominated systems. However, if parts of the reservoir are gravity-dominated either by design or constraint, then density differences in the CO$_\text {2}$-hydrocarbon system become increasingly important. One key effect is the potential for density instabilities to cause gravity override and reduce vertical sweep for systems where CO$_\text {2}$ is lighter than oil (Jenkins 1984), which can be mitigated through alternating CO$_\text {2}$ injection with water, i.e., WAG. There is also potential for convective mixing if density inversions arise in the miscible region (Ahmed et al. 2012; Shahraeeni et al. 2015). Both of these mechanisms may in turn significantly impact the flow behavior of fluids in the mixing zone. While the effect of gravity override is well known, the effect of convective mixing is much less studied.

Density instabilities leading to convection can be created by non-trivial density changes with oil composition. An illustrative example of concentration-dependent density is shown for two binary mixtures, CO$_\text {2}$-butane and CO$_\text {2}$-octane (Fig. 1), as modeled by a cubic equation of state (EOS). First, it is evident that pure CO$_\text {2}$ is lighter than oil at miscible pressures, while a portion of the mixed fluid is heavier than both CO$_\text {2}$ and oil. This results in a non-monotonic density curve shown in Fig. 1a. On the one hand, a light CO$_\text {2}$ fluid will tend to segregate on top of the denser oil, but, on the other hand, diffusion causes CO$_\text {2}$ and oil to mix, creating a dense fluid region in between CO$_\text {2}$ and oil. If this dense mixture evolves on top of the “lighter” oil, an unstable density inversion will occur and lead to density-driven convection. The convective system may be quite complex depending on how the fluids mix. Despite the complexity in density, the viscosity curve (Fig. 1b) remains mostly linear with CO$_\text {2}$ concentration. Significant reduction in oil viscosity will enhance convective mixing, but itself is not an initiating factor in purely gravity-driven setting. We note that for certain high pressure and low temperature reservoirs, CO$_\text {2}$ can be denser than oil, which will lead to different type of density instabilities than described here. Unstable flow occurs in this case when heavy CO$_\text {2}$ is injected above a lighter oil, or from below in combination with other phase-related phenomenon (Moortgat et al. 2013).

Non-monotonicity in density is not unique to CO$_\text {2}$-butane or octane mixtures. Data show that oils exhibit this non-monotonic behavior when mixed with CO$_\text {2}$ under miscible conditions (Lansangan and Smith 1993a, b). However, the exact quantification of these systems in the laboratory remains elusive. The density curves in Fig. 1a have been produced using the Peng–Robinson EOS model. Different thermodynamic models will produce different density curves, i.e., GERG-2008 (Span et al. 2015), and the uncertainty in this regard persists due to the scarcity of data for CO$_\text {2}$-hydrocarbon systems at relevant temperatures and pressures. However, despite the differences between property models, the existing body of evidence suggest that crude oil commonly exhibits similar non-monotonic behavior. Therefore, it is highly likely that non-monotonic density curves are indeed characteristic of CO$_\text {2}$-hydrocarbon systems. Indeed, there are indications of this also at field-scale. For example, in the Weeks Island gravity-stable CO$_\text {2}$ pilot in Louisiana, USA, CO$_\text {2}$ was injected from above a heavier oil in a miscible setting, and breakthrough occurred earlier than expected (Perry 1982). This could have been due to the non-monotonic increase in oil density as it mixes with CO$_\text {2}$, which causes convective mixing. We note that as a general rule, heavier oil components will exhibit greater non-monotonicity. For example, methane will not exhibit such behavior with CO$_\text {2}$ under these conditions, with the density curve for CO$_\text {2}$-methane mixtures being monotonic and close to linear. In addition, the density of lighter oils mixed with CO$_\text {2}$ are more sensitive to pressure and temperature and may become monotonic at lower temperature and pressure.

Convective mixing in a CO$_\text {2}$–brine system is a related and well-studied phenomenon (Elenius et al. 2014; Ennis-King and Paterson 2005; Farajzadeh et al. 2011; Kneafsey and Pruess 2010; Khosrokhavar et al. 2014; Lindeberg and Wessel-Berg 1997; Mykkeltvedt and Nordbotten 2012; Neufeld et al. 2010). In this immiscible two-phase system, density instabilities evolve in the brine phase as CO$_\text {2}$ dissolves across the fluid-fluid interface. Here, CO$_\text {2}$ is assumed to always override the brine due to strong gravity segregation. As CO$_\text {2}$ dissolves from above, brine density increases approximately linearly with CO$_\text {2}$ concentration, causing a density inversion in the brine. Although the density of CO$_\text {2}$-rich brine is heavier than either brine or CO$_\text {2}$, the density differences are slight (up to 1% density increase) and restricted by the solubility limit of CO$_\text {2}$ in brine. Therefore, the density curve for this system is a much simpler monotonic function with respect to CO$_\text {2}$ concentration. There are also no additional destabilizing effects as the impact of dissolved CO$_\text {2}$ on brine viscosity is negligible. Despite small density differences and low solubility, CO$_\text {2}$ dissolution can have a significant impact on large-scale spreading of CO$_\text {2}$ plumes over time scales of 100s or 1000s of years (Gasda et al. 2011, 2012; Elenius et al. 2015).

A number of detailed experimental studies of CO$_\text {2}$–brine convection using analogue fluids have been performed. These systems exhibit non-monotonicity and are of interest for our investigation. For example, mixtures of methanol and ethylene glycol (MEG) with water have a maximum density at intermediate MEG concentration, at certain temperatures and EG fractions (Raad et al. 2016). Raad et al. (2016) show with linear stability analysis (employing the Boussinesq approximation) and direct numerical simulations that the non-monotonicity of this system gives rise to different perturbation growth rates compared to the CO$_\text {2}$–brine system. These analogue systems typically have small density differences (e.g., less than 3% for MEG-water) and account for the lighter “fluid” as a top boundary condition. In contrast, CO$_\text {2}$–oil mixtures can have very large density differences (up to 24% in the systems studied here), may be affected by interaction between the upper and lower region (moving “interface”) that cannot be understood by representation of the lighter fluid as a boundary condition, and may have the lighter fluid initially placed either below or above the heavier fluid, depending on the injection scenario. Previous work by Raad and Hassanzadeh (2015) did, however, employ a moving “interface” between miscible fluids for non-monotonic mixtures of propylene glycol (PPG) and water. This setup was shown to cause slightly faster onset compared to when water was represented by the upper boundary. Non-monotonic density profiles can also develop due to chemical reactions (e.g., Trevelyan et al. 2015) or from mixing of solutes that have composition-dependent diffusivities, or diffusion of both mass and heat (e.g., Trevelyan et al. 2011; Raad et al. 2019). Again, these studies were performed for small density contrasts, with application of the Boussinesq approximation.

In contrast to the systems mentioned above, the potential for density instabilities in CO$_\text {2}$-hydrocarbon mixtures that include large density contrasts and is applicable to a porous media setting has received little attention. Previous studies have focused on field-scale injection of CO$_\text {2}$ into oil using both commercial (Ahmed et al. 2012; Both et al. 2015) and high-resolution simulators (Shahraeeni et al. 2015). This collection of reservoir simulation results demonstrates that gravity fingers develop quickly in both homogeneous and heterogeneous media under realistic reservoir injection conditions. In addition, a few small-scale studies of density-driven convective mixing of CO$_\text {2}$ in oil where there are large density contrasts in different settings: (1) non-monotonic density instabilities in non-porous media (Rongy et al. 2012) and (2) porous media experiments where unstable mixing was due to a viscous drive and other phase partitioning effects (Moortgat et al. 2013). The study of Rongy et al. (2012) confirms that convection occurs due to non-monotonic density effects in their experiments that placed CO$_\text {2}$ above n-decane in a pressure–volume–temperature (PVT) cell, i.e., with no porous medium. Although their gravity-dominated system featured a large density increase of the CO$_\text {2}$–oil mixture and showed that convection accelerated mixing compared to pure diffusive mixing, it was noted by the authors that their analysis of finger growth dynamics was limited to non-porous media. These settings further demonstrate by experiments and simulations that strong convective behavior can be expected where there are large density contrasts. However, additional analysis of the CO$_\text {2}$-hydrocarbon convective system is still needed to understand the impact of non-monotonic density in a porous medium setting.

In our study, we investigate the evolution of gravity instabilities caused by non-monotonicity in oil density as a function of CO$_\text {2}$ concentration. We perform our study using a simple box model to isolate gravity effects from viscous effects. As opposed to previous work, we model not only CO$_\text {2}$ convection from above, but also convection initiated from below. Through detailed numerical simulation, we characterize onset times of convection in both the early (linear) and late (nonlinear) time regimes. Further, we quantify the convective behavior in terms of finger propagation speed and mixing rates, giving insight into the upscaled behavior of convection and its impact on CO$_\text {2}$ uptake in oil. Detailed quantification of these factors has not been performed previously. Importantly, we perform our analysis using a dimensionless simulator that allows our results to be rescaled and applied to other systems of interest. Such generalized analysis is needed to consider field-scale implementations beyond a few specific cases. And finally, the results of our work provide new insights into the potential impact of convection at the field scale.

Fig. 1

Density and viscosity of binary CO$_\text {2}$-hydrocarbon mixtures. Properties calculated for CO$_\text {2}$-butane and for CO$_\text {2}$-octane mixtures at 200 bar and 353 K. We used the Peng–Robinson equation of state for density and linearized viscosity between end points given by Span et al. (2015)

2 Governing Equations and Non-dimensional Formulation

The dynamics of gravitationally driven mixing between fully miscible CO$_\text {2}$ and hydrocarbons in porous media systems can be described in terms of the fluid mass balance, the transport equation, Darcy’s law and constitutional laws for the density and viscosity of the mixture,

$$\begin{aligned}&\phi \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\mathbf {u}}) = 0, \end{aligned}$$

(1)

$$\begin{aligned}&\phi \frac{\partial c}{\partial t} + \nabla \cdot ({\mathbf {u}}c) - \phi D\nabla ^2 c = 0, \end{aligned}$$

(2)

$$\begin{aligned}&{\mathbf {u}}= - \frac{k}{\mu }~(\nabla p -\rho {\mathbf {g}}), \end{aligned}$$

(3)

$$\begin{aligned}&\rho =\rho (p,c), ~ \mu = \mu (c). \end{aligned}$$

(4)

Here, $\rho$ is the mass density [kg fluid/$\hbox {m}^3$ fluid] and c is the molar concentration of CO$_\text {2}$ [mol CO$_\text {2}$/$\hbox {m}^3$ fluid]. In addition, ${\mathbf {u}}$ is the fluid flux [m/s], $\phi$ is the porosity [−], D is the diffusion coefficient [$\hbox {m}^2/\hbox {s}$], p is the fluid pressure [Pa], and ${\mathbf {g}}$ is the gravity vector. It is important to note that we cannot use molar or mass fraction instead of c when modeling concentration-dependent density. In these equations, we have assumed that all species have the same diffusivity and that the mixture is either binary with CO$_\text {2}$ and one hydrocarbon, or the relative composition of hydrocarbon components is constant in time. Therefore, the transport equation is written for CO$_\text {2}$ only. Mechanical dispersion is neglected due to slow velocities of this system. We study a two-dimensional problem with coordinate system (x, z) pointing to the right/up from the initial interface.

To reduce the number of numerical simulations and better understand the governing physics of the mixing process, the system of equations are transformed to a non-dimensional setting. First, we define the following scaled molar concentration and density,

$$\begin{aligned} {{\tilde{c}}}&= \frac{c}{c_\text {max}}, \end{aligned}$$

(5)

$$\begin{aligned} {\tilde{\rho }}&= \frac{\rho -\rho _0}{ \rho _\text {max}-\rho _0} = \frac{\rho -\rho _0}{\varDelta \rho }. \end{aligned}$$

(6)

Here, $c_{{\text {max}}}$ is the molar concentration of pure CO$_\text {2}$ at the pressure at the top of the domain, $\rho _0$ is the density of pure CO$_\text {2}$ at the same conditions and $\rho _{{\text {max}}}$ is the maximum density at this pressure, which is obtained for a mixture of CO$_\text {2}$ and hydrocarbons. For comparison, previous work on dissolved CO$_\text {2}$ in water has often applied a linear density increase with CO$_\text {2}$ concentration (Riaz et al. 2006; Elenius et al. 2014). In that case, the system simplifies such that ${{\tilde{\rho }}}={{\tilde{c}}}$, where the concentration is scaled with its value at the solubility limit. Next, following (Riaz et al. 2006; Elenius et al. 2014) the governing equations are scaled using the following characteristic quantities:

$$\begin{aligned} u^* &= \frac{k_\text {ref}\varDelta \rho g}{\mu _\text {ref}}, \end{aligned}$$

(7)

$$\begin{aligned} l^*&= \frac{D\phi }{u^*}, \end{aligned}$$

(8)

$$\begin{aligned} t^*&= \frac{D\phi ^2}{(u^*)^2}, \end{aligned}$$

(9)

in addition to the following relation between the dimensional and non-dimensional pressure,

$$\begin{aligned} p&= \frac{D\mu _\text {ref}\phi }{k} \left( {{\tilde{p}}}-\frac{\rho _0}{\varDelta \rho } z\right) . \end{aligned}$$

(10)

Here, z refers to the non-dimensional vertical coordinate and $\mu _\text {ref}$ is the reference viscosity, for pure CO$_\text {2}$ at reservoir conditions. In contrast to, e.g., Riaz et al. (2006) and Elenius et al. (2014), the Boussinesq approximation ($\nabla \cdot {\mathbf {u}}=0$) is not applied here due to larger density contrasts. Instead, the full mass balance equation for the fluid is solved. In this case, scaling with the density difference (Eq. (7)) introduces a corrector term $\psi$ in the non-dimensional mass balance equation,

$$\begin{aligned}&\frac{\partial {\tilde{\rho }}}{\partial {\tilde{t}}} + \nabla \cdot (({\tilde{\rho }} +\psi ){{\tilde{{\mathbf {u}}}}}) = 0, \end{aligned}$$

(11)

$$\begin{aligned}&\psi = \frac{\rho _0}{\varDelta \rho }, \end{aligned}$$

(12)

$$\begin{aligned}&\frac{\partial {{\tilde{c}}}}{\partial {{\tilde{t}}}} + \nabla \cdot ({{\tilde{{\mathbf {u}}}}}{{\tilde{c}}}) - \nabla ^2 {{\tilde{c}}}= 0, \end{aligned}$$

(13)

$$\begin{aligned}&{{\tilde{{\mathbf {u}}}}} = - \frac{k/k_\text {ref}}{\mu /\mu _\text {ref}} (\nabla {{\tilde{p}}} + {\tilde{\rho }}), \end{aligned}$$

(14)

$$\begin{aligned}&{\tilde{\rho }}= {\tilde{\rho }}(p,{{\tilde{c}}}), ~\mu =\mu ({{\tilde{c}}}). \end{aligned}$$

(15)

The problem is thus defined in terms of the density and viscosity functions together with initial and boundary conditions. We study homogeneous and isotropic conditions here, i.e., $k/k_\text {ref}=1$.

3 Problem Definition

We study binary mixtures CO$_\text {2}$-C4 (butane) and CO$_\text {2}$-C8 (octane) under isothermal conditions of 353 K. For both mixtures, the dynamics is studied in a square domain with initial placement of a) CO$_\text {2}$ in the upper half (denoted CO$_\text {2}$/C4 and CO$_\text {2}$/C8), and b) CO$_\text {2}$ in the lower half of the domain (denoted C4/CO$_\text {2}$ and C8/CO$_\text {2}$). In each zone, the pressure is initially static, i.e., no flow would occur if this fluid occupied the entire system. We denote the line where the two mixtures initially meet an “interface.” This should not be interpreted as an immiscible interface between phases, because the system is at pressures and temperatures where the components are fully miscible. Further, small perturbations of the concentration with a random distribution between zero and $\varDelta {{\tilde{c}}}=1 \times 10^{-3}$ [−] are initially placed in the grid cells immediately above and below the interface. All boundaries are closed to flow and diffusion except the horizontal boundary of the hydrocarbon component, i.e., the bottom of the domain in placement (a) and the top of the domain in placement (b). This boundary has Dirichlet conditions with initial pressure and zero CO$_\text {2}$ concentration. Thus, flow is allowed through this interface and does occur due to overall density changes.

The characteristic density differences that lead to convection— $\varDelta \rho$, $\varDelta \rho _1$, and $\varDelta \rho _2$—are depicted in Fig. 1 and listed in Table 1. For this initially static problem definition, the driving force behind convection will depend on the initial location of the CO$_\text {2}$. For example, when CO$_\text {2}$ is initially placed above the oil, $\varDelta \rho _2=\rho _{{\text {max}}}-\rho _\text {oil}>0$ (Fig. 1) drives downward propagation of mixed fluid into the initial oil zone. This process depends on diffusion to create the initial mixing zone. In contrast, initial placement of CO$_\text {2}$ below oil leads to two types of convection. First, $\varDelta \rho _1=\rho _{{\text {oil}}}-\rho _{CO_2}>0$ drives upward flow of CO$_\text {2}$ into oil. Second, $\varDelta \rho =\rho _{{\text {max}}}-\rho _{CO_2}>0$ drives downward propagation of dense mixed fluid into the CO$_\text {2}$-rich fluid below.

Table 1 Viscosity and density of pure components and density difference of hydrocarbon-CO$_\text {2}$ mixtures at 200 bars and 353 K

The results are obtained in non-dimensional form but presented in dimensional form for a more direct interpretation. The conversion was made using the scaling factors (7)-(9) presented above, with porosity 0.1, permeability 100 mD ($10^{-13}$ $\hbox {m}^2$), diffusivity $1\times 10^{-9}\,\hbox {m}^2/\hbox {s}$ (Grogan et al. 1988; Renner 1988), in addition to the fluid properties $\varDelta \rho$ and $\mu _{{\text {ref}}}$ given in Table 1.

4 Discretization

We developed a specialized code convect that accurately describes the non-dimensional nonlinear governing equations for simulating convection in porous. A fully implicit formulation solved with a finite-difference discretization is chosen as it is stable, mass conservative, and convergent. Newton’s method is used to linearize the equations, and the Jacobi matrix is determined analytically. We use GMRES with an LU preconditioner to solve the linear sub-problems. The density and viscosity are interpolated from look-up tables, prepared according to the previous section. In the simulations presented here, we use upwinding of advection terms except for the density in Darcy’s law, which is an average to obtain zero flux for static conditions. The non-dimensional pressure is further scaled before solving the system of equations, to reduce the condition number. All simulations employ a 100 x 100 constant-spacing grid.

The current implementation is in MATLAB, where the goal is simulation of fine-scale flow dynamics for simple binary mixtures in small-scale systems as accurately as possible and with the ability to control relevant input parameters and boundary conditions. The purpose is to gain insight into field-scale simulation, such as resolution requirements and possibilities for upscaling, and also to aid in experimental design where boundary conditions are uniquely different from reservoir systems. For simulation of CO$_\text {2}$-hydrocarbon convection in a field setting, any available compositional simulators would be better suited to capture the complexity of multiphase and multicomponent processes such as phase transfer and partitioning. We note that use of higher-order methods has been shown to be preferable for larger-scale simulation of CO$_\text {2}$ convective mixing (Shahraeeni et al. 2015).

The validation of the implementation of the governing equations was performed by several tests, showing that convect could accurately represent the analytical solution to pure diffusion, as well as the dissolution rate of the unstable CO$_\text {2}$–brine system with those density differences and viscosities. The elements of the Jacobi matrix were verified against numerically computed derivatives. For the CO$_\text {2}$–oil systems studied here, the diffusion was again compared to analytical solutions and the small early-time deviations from the diffusion-only profile could be explained by the contraction of the fluids as they mix. Convection was simulated with around 10-20 cells per finger in the lateral direction in most cases. Our earlier studies, e.g., Elenius and Gasda (2013) and Elenius et al. (2014), show that this is sufficient to capture the convection dynamics, and (Elenius et al. 2015) showed that for upscaled parameters such as the dissolution rate, also fewer cells per finger give small errors. Ocular inspection shows reasonable finger shapes. Further validation must await experiments, and the results herein are essential to guide experimental design.

5 Dynamics in the Early Regime

Starting from an initial state of CO$_\text {2}$ above or below the oil with perturbations at the interface, the system evolves due to the combined actions of diffusion and density gradients. Diffusion provides the initial driver for mixing of the fluids, leading to an overall density increase (see Fig. 1). To compensate for this volume reduction, oil is flowing into the domain from the open boundary. As opposed to systems with small density differences such as the CO$_\text {2}$–brine system (i.e., brine with gradients in CO$_\text {2}$ concentration), the inflow of oil initially dominates the flow field. The details of this early regime are described in this section for the CO$_\text {2}$-C8 system based on simulation results.

The simulations presented in this section can be compared with analytical solutions derived from linear stability analysis (LSA). LSA is a common approach employed to determine the critical onset time and wavelengths of the early-time perturbations that lead to convective mixing. In order to perform LSA, the system of equations presented in (11)–(15) must be linearized. This type of analysis has been used successfully to characterize the CO$_\text {2}$–brine convective system as well as non-monotonic systems under the assumption of Boussinesq flow. However, the combination of non-Boussinesq flow together with concentration-dependent viscosity is important for CO$_\text {2}$-hydrocarbon convection, and as such, we were not successful in achieving a solvable system of equations using LSA techniques. Nevertheless, we present the early-time results from simulations to provide data for future LSA studies that are able to overcome this difficulty.

5.1 CO$_\text {2}$ Convection from Above

Figure 2 shows the perturbation of the concentration 1 day and 6 days after initialization of a system with CO$_\text {2}$ placed above the hydrocarbon. The perturbation was calculated by subtracting the horizontally averaged concentration from the concentration. The direction of flow is also indicated. After 1 day, the vertical velocity caused by inflow of oil from below is larger by approximately a factor of 10 than the velocity relating to circular convection at the interface. The dominant finger wavelength close to the interface is a few centimeters, but the pattern is not yet well organized. In the CO$_\text {2}$ region (top), larger eddies form that are shaped by the domain boundaries. After 6 days, the velocity due to inflow of oil is still larger than that of convection at the interface. The perturbations of the concentration have become more organized and now have a defined wavelength of slightly more than 10 cm.

Fig. 2

Streamplot on top of concentration perturbation after a 1 day (non-dimensional time 78000) and b 6 days (non-dimensional time 444000). Non-dimensional times are large because the equations were scaled with $\varDelta \rho$ but the driving force is only $\varDelta \rho _2$

The time chosen in Fig. 2a represents the start of a continuous period of instability in the region below the initial interface ($z<0$). Figure 3 shows the exponential growth rate $\sigma$. As in Elenius et al. (2012), the growth rate is defined based on the rate-of-change of horizontal concentration gradients,

$$\begin{aligned} \sigma&= \frac{1}{\varDelta t}\ln \left( \frac{\int | \partial c/\partial x |(t) dx dy}{\int | \partial c/\partial x |(t-\varDelta t) dx dy}\right) . \end{aligned}$$

(16)

The regions above and below the interface are computed separately into two values of $\sigma$. Instability begins in the lower region where the density profile is internally unstable, increasing up toward the interface. There is an early period of varying stability and instability above the interface. The inflow of oil due to density changes around the interface also has an impact on the perturbation growth. Here, the linear onset time $t_c$ in either of the two regions is defined as the time from which the system is continuously unstable in that region. We have $t_c\approx 1$ day below the interface and $t_c\approx 2$ days above the interface. Note that also the perturbation growth above the interface is due to the unstable profile below the interface by continuity of flow patterns. The onset of instability 1–2 days is very quick, and for any practical purpose it can be said to be instantaneous. It should be noted that interpretation of simulation results are based on all nodes in an average sense, as opposed to the most unstable mode used in LSA.

Fig. 3

Exponential growth rate $\sigma$ when CO$_\text {2}$ is placed above the oil. The lower and upper regions are continuously unstable from $t\approx 1$ day and 2 days, respectively

If there were no density differences, the early diffusive profile would be given by

$$\begin{aligned} c_0(\xi ,t)&= {\bar{c}} - (\varDelta c/2)~\mathrm {erf}(\xi ). \end{aligned}$$

(17)

Here, $\xi =z/(2\sqrt{t})$ where z is the non-dimensional distance above the interface and t is the non-dimensional time. We also define ${\bar{c}}=(c_l+c_u)/2$ and $\varDelta c=c_l-c_u$ with $c_l$ and $c_u$ representing the initial concentrations in the lower and upper regions, respectively (here having the values 0 or 1). This diffusion profile is shown at 1 day by a dashed line in Fig. 4a. The simulation result (black solid line) deviates slightly from the pure-diffusion profile due to density changes. Note that the vertical position is given in centimeters above and below the interface on the left axis and in terms of the self-similar coordinate for pure diffusion, $\xi$, on the right axis. Figure 4b shows the corresponding results for the density profile in black. The density profile with maximum density close to scaled concentration 0.5 (Fig. 1) is consistent with a density maximum at the initial interface with an internally stable density profile above and an internally unstable density profile below. This is consistent with instabilities starting in the lower region, as discussed above. The largest perturbations have not yet shifted from the initial position at the interface although the density gradient is more adverse further down. Figure 4 (right panels) also shows corresponding results at 6 days. The diffusion profile is wider, and the position of maximum growth has shifted downward toward a more unstable density gradient, consistent with what is shown in previous figures.

Fig. 4

Simulation results after 1 day (left panels) and 6 days (right panels). Top: horizontally averaged concentration. Bottom: horizontally averaged density and nominal perturbation size

5.2 CO$_\text {2}$ Convection from Below

Figure 5 shows the perturbation of the concentration and the flow direction at 122 s and 731 s, respectively, after initialization of a system with CO$_\text {2}$ placed below the hydrocarbon. After 122 s, the vertical velocity caused by inflow of oil from above is approximately a factor 40 larger than the velocity relating to circular convection at the interface. The dominant finger wavelength close to the interface is approximately 2 mm, corresponding to a non-dimensional length around 50, which is 5 grid cells. Further down, larger eddies form that are shaped by the domain boundaries. There is flow toward the lower left corner but with velocities tending to zero.

After 731 s, the velocity due to inflow of oil is still larger than that of convection at the interface, but now it differs only by a factor of two. The perturbations of the concentration have the same horizontal wavelength as before but the largest perturbations occur approximately 2 mm below the interface. The flow field below the fingers has become more organized with three modes, except at the bottom of the domain where there is only one mode.

Fig. 5

Streamplot on top of concentration perturbation after a 122 s (non-dimensional time 100) and b 731 s (non-dimensional time 600)

The time chosen in Figure 5a represents the start of a continuous period of instability. Figure 6 shows the exponential growth rate $\sigma$. As before, the regions above and below the interface are computed separately into two values of $\sigma$. The upper region is stable at early times. Instead, instability begins in the lower region where the density profile is internally unstable, increasing up toward the interface. Similar to the case with CO$_\text {2}$ above the hydrocarbon, the full system features an early period of varying stability and instability. The inflow of oil due to density changes also has an impact on the perturbation growth. We observe an onset time of $t_c\approx 100$ s below the interface and $t_c\approx 200$ s above the interface. The onset of instability is therefore extremely quick and for any practical purpose it can be said to be instantaneous.

Fig. 6

Exponential growth rate $\sigma$ when CO$_\text {2}$ is placed below the oil. The lower and upper regions are continuously unstable from $t\approx 100\,\hbox {s}$, respectively, 200 s

The concentration profile due to pure diffusion is shown at 122 s by a dashed line in Fig. 7a. The simulation results (black solid line) deviate slightly from the pure-diffusion profiles due to density changes. Figure 7b shows the corresponding results for the density profile in black. The density has a maximum at the initial interface with an internally stable density profile above and an internally unstable density profile below. This is consistent with instabilities starting in the lower region, as discussed above. The largest perturbations (red line) occur less than 1 mm below the interface at this time, where the density gradient is most adverse. Figure 7 (right panels) also shows corresponding results at 731 s. The diffusion profile is wider, and the position of maximum growth has shifted downward, consistent with what is shown in previous figures. The diffusion-only profile gives a worse estimate at this time.

Fig. 7

Simulation results after 122 s (left panels) and 731 s (right panels) Top: horizontally averaged concentration. Bottom: horizontally averaged density and nominal perturbation size

In summary, an initial placement of CO$_\text {2}$ either above or below octane is very unstable with onset occurring within a few days. The most unstable setting is placement of CO$_\text {2}$ below the hydrocarbon.

6 Dynamics in the Fully Nonlinear Regime

6.1 CO$_\text {2}$ Convection from Above

We begin by examining the situation with initial CO$_\text {2}$ placement above the hydrocarbon. Figure 8 shows the scaled CO$_\text {2}$ concentration and the fluid density resulting from CO$_\text {2}$ mixing with C4 (left sub-figures) and with C8 (right sub-figures). Note that the domain scaled to the chosen parameters is less than a meter in each direction, and that the finger width is on the scale of 10 cm. Thus, it is necessary to employ sub-centimeter grid resolution to simulate the fingers accurately. As the fingers form and migrate downward, the initial interface between the CO$_\text {2}$ and hydrocarbon retreats in time due to mixing of CO$_\text {2}$ to the oil region. The arrows show the direction of flow. There is inflow of oil across the open lower boundary. This is compensating for the decrease in fluid volume as the components mix (increased density). Just below fingers, there is also outflow of oil across the boundary, following the flow direction of fingers.

Fig. 8

CO$_\text {2}$ above C4 (left) and C8 (right). For each system, we show the scaled CO$_\text {2}$ concentration [-] and the density [kg/$\hbox {m}^3$]

The fingers do not form immediately. As in Elenius and Johannsen (2012), we define the nonlinear onset time $t_{{\text {nonlin}}}$ as the time when fingers can first be observed (without subtracting the horizontally averaged concentration) and when convection has a visible impact on the rate of transfer of components between zones. Table 2 presents this onset time together with the speed of fingers and the speed of the retreating interface, while Fig. 9a provides the evolution of the fingers and retreating interface in time. The onset time is approximately 3 days for CO$_\text {2}$ convective mixing in C4 and 22 days for convective mixing into C8. The latter can be compared with the linear onset time (start of perturbation growth) which is 1 day. The downward finger speed after onset is 4 cm/day for C4 and only 2 cm/day for C8, despite very similar $\varDelta \rho _2$. This is explained by differences in viscosity. Although the CO$_\text {2}$ viscosity is the same in both systems, the viscosity of C8 is a factor 3 larger than the viscosity of C4.

The lower sub-figure of Fig. 9 shows the mass transfer rate of CO$_\text {2}$ in time for both mixtures varies steadily over time until the fingers reach the bottom of the domain. The system evolves quickly, but before the fingers reach the bottom, the long-term average is around 0.3 kg/$\hbox {m}^2$day for C4 and 0.1 kg/$\hbox {m}^2$day for C8. For comparison, the dissolution rate of CO$_\text {2}$ into brine with a sharp-interface assumption is $5 \times 10^{-4}$ kg/$\hbox {m}^2$day at a similar permeability and salinity from CO$_\text {2}$ injection into Tubåen at the Snøhvit CCS project, which has similar reservoir conditions. The transfer rates for C4 and C8 are more than 100 times larger than that of CO$_\text {2}$ mass transfer into brine. This large increase is due to a combination of increased density difference and reduced viscosity in a fully miscible setting for CO$_\text {2}$-hydrocarbon versus CO$_\text {2}$–brine mixtures.

Table 2 Onset time of convection (i.e., nonlinear onset time), finger speed and mixing rate in systems with CO$_\text {2}$ initially placed above, respectively, below, the hydrocarbon

Fig. 9

CO$_\text {2}$ above hydrocarbon. Above) Advancing finger tips defined by lowermost occurrence of ${{\tilde{c}}}>0.2$ for CO$_\text {2}$ fingers migrating down (red) and uppermost occurrence of ${{\tilde{c}}}<0.8$ for hydrocarbon front migrating up (black). Blue and burgundy lines mark the location of domain boundaries. Due to the open bottom boundary, fingers leave the domain after 13 days (C4) and 34 days (C8). Below) Uptake rate of CO$_\text {2}$ into oil region below (defined as having ${{\tilde{c}}}<0.8$)

6.2 CO$_\text {2}$ Convection from Below

We now examine the situation when CO$_\text {2}$ is initially placed below the hydrocarbon. Because this system differs more from previously studied CO$_\text {2}$–brine systems, we will examine the dynamics in more detail. Figure 10 shows the scaled CO$_\text {2}$ concentration at different times up to 16 days, when CO$_\text {2}$ is placed below C4. The corresponding density of the fluid is also shown. Again, when rescaling to the chosen parameter set, the finger widths are on the centimeter scale. The initial system is almost immediately unstable since CO$_\text {2}$ is lighter than the oil, and we see fingers of CO$_\text {2}$ migrating up into the oil. Because of diffusion, intermediate concentrations are obtained along the outer edge of the fingers, shown in yellow-green shading, which have a higher density than either the oil or the CO$_\text {2}$. As a consequence, the mixed fluid flows down into the CO$_\text {2}$ region, also in the shape of fingers.

The downward migration of mixed fluid into the CO$_\text {2}$ is much faster than upward migration of CO$_\text {2}$. The reason for this is that the driving force of downward migration is $\varDelta \rho >\varDelta \rho _1$, cf. Fig. 1 and Table 1. In addition, the CO$_\text {2}$ fingers migrating upward are being consumed as they propagate, which does not happen for the fingers moving down into the CO$_\text {2}$ region.

Since the initial CO$_\text {2}$ layer has a finite thickness, mixing in the lower region eventually reduces the CO$_\text {2}$ concentration there and thus reduces its buoyant force relative to the upper region. In this example with a CO$_\text {2}$ layer of 42 cm, the lower region is no longer buoyant with respect to the oil after 8 days. After this time, there is only internal convection within the lower region until the concentration there is homogeneous. Perhaps, the feature of largest practical importance in this system is that there is no more CO$_\text {2}$ migration up into the oil by convection after this time. The figure also shows that oil enters the domain through the open top boundary, again due to the overall density increase as the fluids mix.

Fig. 10

Scaled CO$_\text {2}$ concentration (top) and density (bottom) when CO$_\text {2}$ is placed below C4

Fig. 11

Scaled CO$_\text {2}$ concentration (top) and density (bottom) when CO$_\text {2}$ is placed below C8

Figure 11 shows the same dynamics evolve when CO$_\text {2}$ is placed below C8. Table 2 presents the speed of fingers and the rate of oil transfer into the CO$_\text {2}$. This is shown in more detail in Fig. 12. For all practical purposes, the migration of fingers when CO$_\text {2}$ is placed below the oil can be said to begin immediately after initialization of the system. The left sub-figure shows the location of the lower-most and upper-most finger tips. The speed of downward migrating fingers is 20 cm/day in both systems, and that of upward migrating CO$_\text {2}$ fingers is 5 cm/year into C4 and 10 cm/year into C8. The increased speed of finger migration into C8 (despite being more viscous) is caused by the considerably larger $\varDelta \rho _1$ in this system (Fig. 1 and Table 1).

The right sub-figure of Fig. 12 shows the mass transfer rate of hydrocarbon into the CO$_\text {2}$ plume. This is an order of magnitude larger than the mass transfer of CO$_\text {2}$ when convecting from above. The mixing of oil into CO$_\text {2}$ is therefore very efficient when CO$_\text {2}$ originates from below.

Fig. 12

CO$_\text {2}$ below hydrocarbon. Above) Advancing finger tips defined by uppermost occurrence of ${{\tilde{c}}}>0.8$ for CO$_\text {2}$ fingers migrating up (red) and lowermost occurrence of ${{\tilde{c}}}<0.2$ for hydrocarbon fingers migrating down (black). Blue and burgundy lines mark the location of domain boundaries. Below) Uptake rate of oil in CO$_\text {2}$ region below (defined as having $c>0.2$)

7 Field-scale Implications

We have studied systems that are purely gravity driven within a simple box model. Although our study is not intended to exactly replicate complex flowing conditions of EOR or CO$_\text {2}$ storage operations, we can nevertheless use our results to gain important insights into field-scale impacts of convective mixing during a miscible flood for systems where CO$_\text {2}$ is lighter than oil. Our simple system may be considered an approximation of regions in the reservoir where gravity is more dominant than viscous forces. This may be the case in reservoirs with large distances between wells, or in the case of CO$_\text {2}$ injection into depleted oil and gas reservoirs where mixing of CO$_\text {2}$ and remaining hydrocarbons can occur.

One good example is the Snøhvit CCS project where CO$_\text {2}$ has been injected for several years into the Stø formation that is an active gas-producing reservoir (Kaufmann and Skurtveit 2018). Although this case is a storage project and not for EOR, CO$_\text {2}$ was injected in the water leg in the vicinity of a 10-m oil zones, and thus, there are likely interactions between the injected CO$_\text {2}$ and resident oil and gas. In fact, there could be risk of the gas production being contaminated by CO$_\text {2}$. Our study indicates that CO$_\text {2}$ breakthrough into the gas cap would have to migrate upward through the oil. We show that the buoyancy drive of CO$_\text {2}$ would be significantly delayed by downward convection of a dense CO$_\text {2}$–oil mixture. If the system in the reservoir were purely gravity driven, we estimate that the accumulation of CO$_\text {2}$ below the oil would need to be 30-m thick before CO$_\text {2}$ could counter the downward convection and emerge into the gas cap. We caution that since we are not using the exact oil properties of the Stø reservoir or simulating the full 3-phase system, our estimate is very approximate. However, the results of our studies implies that CO$_\text {2}$ injection into Stø needs to be able to capture the potentially important convective processes.

Another example where convection may play a role is in cases of significant gravity override during a CO$_\text {2}$-flood. Gravity override is understood as reducing reservoir sweep and early CO$_\text {2}$ breakthrough. Here, convection acts to dissolve the CO$_\text {2}$ from above into the oil below. If a 1-m-thick gravity tongue of CO$_\text {2}$-rich fluid evolves over several meters of oil, then our results suggest that convection would reduce the plume thickness at a rate of 0.3–1 cm/day, resulting in completely dissolved CO$_\text {2}$ within 3 months to 1 year. We emphasize that the exact dynamics of CO$_\text {2}$ dissolution needs to be investigated under flowing conditions to understand the interplay between viscous and gravity effects. Characterization of the competition between viscous- and gravity fingering is currently an active area of research, but it implies that a gravity tongue would likely migrate more slowly through the reservoir under the influence of convection. In addition, reservoir sweep could be enhanced as CO$_\text {2}$ convects through oil to deeper parts of the reservoir.

8 Discussion

We have presented results and analysis of convection in CO$_\text {2}$-hydrocarbon systems using highly resolved simulation. The system presented herein is a relatively simple box model with a two-component single-phase fluid under initially static conditions but where evolution of CO$_\text {2}$-hydrocarbon mixing is surprisingly complex. The driving forces are density instabilities arising from the non-monotonic nature of oil density with increasing CO$_\text {2}$ concentration that exhibits a maximum mixture density occurring at intermediate CO$_\text {2}$ mole fraction. Viscosity reduction with CO$_\text {2}$ dissolution further enhances the convection initiated by density instabilities.

An important finding of this work is characterization of temporal and spatial scales of the convective process, from early- to late-time regimes. The simulator is implemented in non-dimensional form, meaning that length and time scales of our results can be easily rescaled to suit the reservoir of interest by applying Eqs. (7)–(9) with the desired parameters. In this study, we have chosen conditions typical of oil and CO$_\text {2}$ storage reservoirs (100 mD). One may need to re-simulate convection if the corrector term in Eq. (12) changes significantly, i.e., a different oil composition. For our chosen parameter set, the width of convective fingers is on the order of centimeters. However, an increase in permeability by a factor of 10 (to 1 Darcy) would result in a factor of 10 decrease in finger width (to millimeters), and a factor of 100 decrease in the time scale (from days to minutes).

Consistent with previous work, our results show that accurate simulation of convective mixing for typical reservoirs requires very high-resolution grids, with millimeter-sized cells needed for centimeter-scale fingers. This implies that direct simulation of convective mixing is not possible at the reservoir scale ($100\, \hbox {km}^2$) using any type of compositional simulator. We note that if fingers are larger (meter-scale) (Ahmed et al. 2012), then it may be worthwhile employing high-resolution techniques following previous work of Shahraeeni et al. (2015). In general, new simulation techniques are needed to model the impact of sub-meter-scale convection at the scale of typical field-scale cell sizes, and the resolved box-model simulations presented herein can guide development of upscaled models. Another alternative would be to employ dynamic grid refinement approaches. Both subjects are of great interest, but currently beyond the scope of this work.

Experiments of CO$_\text {2}$-oleic mixing in a controlled setting are needed to verify the findings of this paper. Our results will be helpful to guide experimental setup, giving expected flow behavior and providing scaling relations to find the correct length and time scales of the experiment. We recommend that the dimensions and flow properties of the experimental setup be chosen to include at least a few fingers in the lateral direction, and imaging should provide a resolution of a several pixels per finger. In addition, the height of the domain should be large enough to be able to observe the regime before fingers interact with the bottom. The setup must allow for inflow of oil to compensate for the pressure reduction as the fluids mix. Preferably, the inflow boundary should be situated at a distance from the initial interface to reduce perturbations. Since reliable thermodynamic and other fluid data are still lacking for CO$_\text {2}$-hydrocarbon mixtures, measurements of density, viscosity and diffusivity would be highly useful for calibrating models.

The study presented herein presents important insights into the implications of fine-scale convection in CO$_\text {2}$-hydrocarbon systems, which is a relatively less studied effect than viscous-capillary effects in CO$_\text {2}$ flooding or compared to convection in CO$_\text {2}$–brine systems. As such, we have not addressed all theoretical, numerical or practical aspects of convection for these systems that would be interesting areas for future work. We hope that a rich body of literature can be developed for CO$_\text {2}$–oil convection as has occurred for CO$_\text {2}$–brine systems. For example, we have focused on pressure and temperatures relevant for CO$_\text {2}$-EOR to obtain the density and viscosity curves for the fluid pairings used in this study, but it would be interesting to explore the impact of density separate from viscosity for abstract fluid pairings. Another area of study could be to investigate the impact of the shape of the density curve rather than only the extreme points (as shown by (Raad et al. 2016) in the context of MEG solutions), which could be addressed by including initial and boundary conditions to feature density maximum at smaller and larger CO$_\text {2}$ concentrations than included here. And finally, one could investigate the impact of concentration-dependent diffusivity, a known effect in hydrocarbon systems (Leahy-Dios and Firoozabadi 2007), which would necessitate a fully compositional simulator and not the relatively simple simulator developed for this study. In conclusion, there are a range of aspects to CO$_\text {2}$-hydrocarbon convection that can and should be further explored in future work.

9 Summary and Conclusions

Our study has demonstrated that convective mixing of CO$_\text {2}$ into oil in our examples is highly efficient. The onset time and rate of mixing depends on the position of the CO$_\text {2}$-rich fluid compared to the oil. Convection from above is caused by an increase in oil density at low CO$_\text {2}$ concentrations that initiates a density inversion in the oil and leads to convection. Conversely, convection from below is a more complex phenomenon where two-way convection occurs as CO$_\text {2}$-rich fluid rises and is countered by sinking mixed fluid. We find that under typical reservoir conditions, convection starts within few days if CO$_\text {2}$ is introduced from above and almost immediately when CO$_\text {2}$ convects from below. CO$_\text {2}$ propagates through the oil quickly, reaching a meter within 1–2 weeks in the case when convection originates from below, and up to a month when originating from above.

The onset time and convection rates for the CO$_\text {2}$-hydrocarbon systems studied here are far more efficient than those of the well-studied CO$_\text {2}$–brine system under similar conditions. We find convection rates on the order of 200–600 times greater than CO$_\text {2}$ convection in brine. This result suggests that convective mixing in CO$_\text {2}$-hydrocarbon systems is a process that occurs at relevant timescales of CO$_\text {2}$ injection, i.e., days to months, compared to the CO$_\text {2}$–brine system where convection acts over much longer timescales, i.e., decades to centuries. Therefore, an important conclusion of this work is that the possibility for convection should be considered in CO$_\text {2}$-EOR or CO$_\text {2}$ storage projects where hydrocarbons are present.