Numerical investigation of convective mixing in impure CO2 geological storage into deep saline aquifers

International Journal of Greenhouse Gas Control 96 (2020) 103015

Contents lists available at ScienceDirect

International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc

Numerical investigation of convective mixing in impure CO2 geological storage into deep saline aquifers

T

Didi Lia, Xi Jiangb,* a

Guangzhou Key Laboratory of Environmental Catalysis and Pollution Control, School of Environmental Science and Engineering, Institute of Environmental Health and Pollution Control, Guangdong University of Technology, Guangzhou, 510006, China b School of Engineering and Materials Science, Queen Mary University of London, Mile End Road, London, E1 4NS, UK

A R T I C LE I N FO

A B S T R A C T

Keywords: CO2 geological storage Impurity Convective mixing Scaling relation Numerical simulation

Dissolution of CO2 into formation brine could lead to density-driven convective mixing which enhances the transfer rate of free-phase CO2 into saline aquifers and thus plays an important role in efficient and secure longterm geological storage. In practical projects, CO2 is always co-injected with multiple impurities to reduce capture costs or to dispose hazardous species like SO2 and H2S alongside. Numerical simulations in this study reveal that, compared with the pure CO2 case, only the co-injection of SO2 would strengthen the convective mixing process in the formation brine while all the other common impurities would result in delayed onset and weakened convection. Particularly, single impurity such as N2 may represent all the air-derived impurities in predicting the evolution of convective mixing process. With the exception of SO2, the diffusivity of most coinjected impurities is not likely to have noticeable impact on the occurrence and development of the densitydriven convection. However, the effective CO2 diffusivity should be measured or estimated when SO2 is coinjected. Several scaling relations characterizing the onset, strength, and decay of the convective mixing process in impure CO2 geological storage are also established, which would provide reference for screening potential reservoirs for impure CO2 geological storage.

1. Introduction CO2 geological storage (CGS) into deep saline aquifers is considered as one of the most promising options to reduce net CO2 emissions into the atmosphere. To ensure the effectiveness of CGS and finally meet the fundamental aim of mitigating global warming, the retention time of CO2 in the subsurface is supposed to be in the order of several thousand years or more with annual leakage rates not exceeding 0.01% (Alcalde et al., 2018; Metz et al., 2005). It is thus important to fully understand the trapping mechanisms in the saline aquifers to minimize the risk of leakage. There are mainly four trapping mechanisms in the saline reservoirs, including static trapping, residual-gas trapping, dissolution trapping and mineral trapping (Bachu, 2008; Chadwick et al., 2008; Jiang, 2011). CO2 dissolution trapping occurs when injected CO2 dissolves and becomes trapped within the formation brine. This mechanism plays an important role in the security of CGS in the medium to long term. More importantly, CO2 dissolution is the prerequisite of the permanent and more secure mineral trapping that occurs over the 1000-year timescale (Alcalde et al., 2018). After injected into the reservoir, the CO2 stream tends to rise

⁎

upwards due to buoyancy forces and accumulate under the caprock, forming a layer of free-phase CO2 and starting to diffuse into the underlying aqueous phase, creating a diffusive boundary layer which grows according to the Fickian diffusion law (Emami-Meybodi et al., 2015). Depending on pressure, temperature and salinity of the formation brine, dissolved CO2 increases the density of the aqueous phase by about 0.1%–1%, generating a gravitational instability, which could give rise to density-driven convective mixing which accelerates the dissolution rate of free-phase CO2 into the formation brine (EmamiMeybodi et al., 2015; Taheri et al., 2018). It is necessary to make precise prediction of the onset time and the strength of this convective mixing because the timescale for dissolution corresponds to the timescale over which free-phase CO2 has a chance to leak out from the storage sites. It should be noted that although the injected stream is “overwhelmingly CO2” (Gale, 2009), certain impurities such as N2, O2, Ar, H2S, SO2, and CH4 may also be co-injected. Some impurities are allowed to be included in order to reduce the cost of separation and compression during the capture process while hazardous ones like H2S or SO2 may be intentionally co-injected to be disposed of (Crandell

Corresponding author. E-mail address: [email protected] (X. Jiang).

https://doi.org/10.1016/j.ijggc.2020.103015 Received 29 September 2019; Received in revised form 9 February 2020; Accepted 5 March 2020 Available online 13 March 2020 1750-5836/ © 2020 Elsevier Ltd. All rights reserved.

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

Nomenclature

z

c0 c1 D g g H k K l L M n P q Ra Sh t Vi Vn W x y

Greek symbols

numerical constant for onset time numerical constant for critical wavelength gas diffusivity in the aqueous phase (m2/s) magnitude of gravitational acceleration (m/s2) vector of gravitational acceleration (m/s2) height of the simulation domain (m) formation permeability (m2) equilibrium constant thickness of the diffusive boundary layer (m) length of the simulation domain (m) molecular weight normal vector pressure (bar) Darcy velocity (m/s) Rayleigh number Sherwood number time (s) partial molar volume of the dissolved gas (m3/mol) mesh volume (m3) width of the simulation domain (m) component mole fraction in the aqueous phase component mole fraction in the gas phase

β ϕ λ μ ρ ρo τ

component mole fraction in the total system

mole fraction of gas phase in the system porosity wavelength (m) dynamic viscosity (Pa∙s) molar density (mol/m3) mass density (kg/m3) dimensionless time

Subscripts/superscripts aq c d i imp m max n w

aqueous phase critical decay component index impurity number of components maximum grid block index water

develop slower mixing rate (Raad and Hassanzadeh, 2017). Meanwhile, using both linear and nonlinear analyses, Kim and Song (2017) implied that double diffusive effects may have positive or negative effects on the system stability based on the values of the diffusivity ratio and buoyancy ratio, e.g., impurities such as H2S and N2 may have the potential to retard the onset of instability and reduce the dissolution rate. The critical difference between these two sets of numerical results may arise from the different definitions of the effective Rayleigh number (Raad and Hassanzadeh, 2018; Kim and Song, 2019). Mahmoodour et al. (2018) estimated the effective CO2 diffusion coefficient in the CO2-N2brine system to calculate the effective Rayleigh number based on total flux. Their experimental and theoretical results suggested that the inclusion of the N2 impurity may accelerate or delay the onset of convection depending on the concentration of N2. For example, the 10% N2 would promote the onset time while the 20% N2 resulted in a similar onset as that of the pure CO2 case. The aforementioned experiments, simulations, and theoretical analyses contribute to a better understanding of the convective mixing process in impure CO2 storage. However, many studies were focused on only one specific impurity. Even for one specific impurity, there are contradictory conclusions resulting from different research methods and definitions of characteristic numbers. Furthermore, studies on the effects of the diffusion coefficients of the co-injected impurities on the convective mixing are limited. It seems necessary to obtain the effective CO2 diffusion coefficient in order to calculate the effective Rayleigh number for impure CO2 geological storage. However, direct measurements of the gas diffusivity in the formation fluids are rather complicated and require special equipment while indirect estimations may lead to uncertainties. The present study aims to analyze and compare the effects of different common impurities, including N2, Ar, CH4, H2S, O2 and SO2. We would also like to investigate the potential of single impurity as the representative for multiple impurities in the convective mixing studies. In addition, we intend to examine the effects of impurity diffusivity on the onset time and development of the densitydriven convection process, aiming at determining the necessities of the measurements or calculation of the effective CO2 diffusion coefficient in certain impure CO2-solution system. More importantly, scaling relations to predict the evolution of the convective mixing process in impure CO2 storage would be obtained to help select suitable storage sites,

et al., 2010; Jacquemet et al., 2007). The presence of non−CO2 species is supposed to affect the properties of the CO2 stream and thus the trapping mechanisms in the reservoir (Talman, 2015). In the first place, non-condensable impurities, especially light-impurity species such as N2 and CH4, would increase the buoyancy force for the CO2 plume, reducing the efficiency of CO2 dissolution into the formation brine (Wang et al., 2011a, b). Apart from the decreased density of the CO2 stream, most impurities were implied to reduce the density of the aqueous phase when dissolved at reservoir conditions (Ennis-King and Paterson, 2005). For example, H2S was suggested to have negative effects on the density of the aqueous solution. It is thus predicted that the inclusion of H2S was unfavorable for the density-driven convection and thus the dissolution trapping mechanism (Ji and Zhu, 2013). Further numerical simulations of the convective mixing process confirmed the inhibitive effects of H2S and provided insights into the optimal H2S composition in the CO2-H2S mixture in geological disposal (Liu et al., 2017). Furthermore, N2 impurity was suggested to result in delayed onset time of convection and reduced dissolution rate while SO2 impurity was implied to enhance the density-driven convective mixing process (Li and Jiang, 2014; Li et al., 2015; Mahmoodpour et al., 2017). However, these investigations assumed equal diffusivities of CO2 and impurities, ignoring the double diffusive effects (Kim and Song, 2017; Raad and Hassanzadeh, 2016). In fact, the presence of impurities may change the diffusion coefficient of CO2. For instance, diffusion coefficient of CO2 in the CO2-N2 mixture was supposed to be smaller than that of pure CO2 alone (Mahmoodpour et al., 2018). Taking into consideration the effects of different diffusivities of CO2 and impurities, it is suggested that the H2S or N2 impurity has both the possibilities to impede and accelerate the density-driven convection, depending on the composition of the injected impure CO2 stream (Kim and Song, 2017, 2019; Mahmoodpour et al., 2018; Raad and Hassanzadeh, 2016, 2017, 2018). To be more specific, Raad and Hassanzadeh (2016) conducted linear stability analyses and suggested that the diffusion contrast between CO2 and the impurity species such as H2S would result in a nonmonotonic density profile which may lead to an earlier onset of convective instability and a possible enhancement of the solubility trapping mechanism. Furthermore, their simulation results indicated that the convective mixing dynamics of an impure CO2-H2S system may change dramatically after the onset, i.e., scenarios with earlier onset may 2

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

of all gases are set as 2 × 10−9 m2/s. In Section 4.3, we further explore the effects of diffusion coefficients of different gases. The governing equations used in the present study are derived from conservation laws (Li and Jiang, 2014; Li et al., 2015). Since nonthermal condition is considered in the present study, the energy conservation equation is neglected. Besides, there is no source or sink term in the model. As a result, the mass conservation equations specified to our multi-component aqueous phase model can be written as follows,

which can be of practical significance. 2. Model setup In a typical geological saline aquifer, the buoyant CO2 mixture injected into the permeable storage reservoir would displace the formation fluids and accumulate beneath the impermeable caprock. A threedimensional domain below the two-phase interface is chosen as the research domain (Fig. 1). It is worth mentioning several relevant features and assumptions of the present model. First, the two-phase interface is supposed to remain relatively sharp and will not be deflected by the convective mixing in the underlying aqueous phase (Farajzadeh et al., 2007; Li and Jiang, 2014). Second, since this study is mainly focused on the fluid flow and transfer in the aqueous phase and instead of dealing with complicated two-phase flows involving potential phase transition as well as capillary pressure and relative permeability, it is assumed that the model is comprised of single aqueous phase and the top boundary of the model is represented as gas-saturated water. For the length of the simulation, the free gas mixtures above the top boundary are supposed to be abundant to keep the aqueous phase at the top boundary saturated with dissolved gases. Third, it has been suggested that the integral measures, including the onset time of convection and the dissolution rate, are not sensitive to the variations of the lateral boundary conditions as well as the bottom boundary conditions (Zhang et al., 2011). As a consequence, “no flow” boundary conditions are adopted on these boundaries for simplicity. Fourth, the initial ambient fluid in the simulation domain is assumed to be pure water, with no gases dissolved. Fifth, temperature is assumed to be constant for the small-scale model. Last but not least, in real-world geologic reservoirs, Rayleigh instabilities are often triggered by geologic heterogeneities such as heterogeneities in permeability distribution (Xie et al., 2011). In the present study, a small random fluctuation (up to ± 4%) is imposed on the permeability field to mimic the real reservoirs as well as to trigger the onset of convection. The parameters used in the present investigations, including the fluid and formation properties as well as the model size, are specified in Table 1. Six kinds of common non−CO2 species are considered, including Ar, CH4, H2S, N2, O2 and SO2. The upper limit of impurity mole fraction in the injected CO2 mixture is 10% because it is the possible maximum of the impurity concentration in the captured streams (Li et al., 2011). In order to compare and analyze the effects of different impurities at the same concentrations, the impurity concentrations in the initial feed gas compositions are all chosen to be the same as 10%. By assuming local thermodynamic equilibrium, vapor-aqueous flash calculation is then used to calculate the saturated gas compositions in the aqueous phase at the top boundary from the injected feed gas compositions (Li and Jiang, 2014; Ziabakhsh-Ganji and Kooi, 2012). In this work, the Rachford-Rice equation is used to obtain the partitioning of different components into two phases (Rachford and Rice, 1952; Leibovici & Neoschil, 1995), given as

z (K − 1) =0 i − 1)

∑ 1 +i β (iK i

∂ ∂t

∫ ϕρxi dVn + ∫ ∇⋅(ρxi q) dVn = ∫ ∇⋅(ϕDρ∇xi ) dVn Vn

Vn

Vn

k q = − (∇P − ρo g ) μ

(3)

(4)

with the supplementary constraints as, m

∑ xi = 1

(5)

i=1

Maq ρaq

=

x w Mw + ρw

m−1

∑ Vi xi i=1

(6)

Eq. (6) is used to calculate the aqueous phase density (ρaq) by taking into consideration the effects of dissolved gases, where m-1 stands for the number of components excluding water. The partial molar volume of dissolved CO2 (VCO2) is dependent on temperature only (Garcia, 2001) while the partial molar volume of dissolved N2 is related to both temperature and pressure (Mao and Duan, 2006). The SUPCRT92 package, a software package for calculating the thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0–1000 °C (Johnson et al., 1992), is adopted to compute the partial molar volume of other investigated gases. No attempt is made to take into consideration the effects of dissolved gases on the viscosity of the aqueous phase, i.e., the viscosity of the aqueous phase is constant throughout the simulation. In order to capture the early-time diffusive behavior, grid sizes near the top boundary should be smaller than the critical thickness of the diffusive boundary layer at the onset of convection, ϕμD lc = 3.6 c0 k Δ ρg = 0.011m (Ennis-King et al., 2005). Away from the top boundary, the vertical grid spacing can be gradually increased in order to improve computational efficiency. For the other two dimensions, the grid sizes are supposed to be at least sufficient to resolve the critical finger wavelength at the initiation of the convective instability ϕμD (λ c = c1 k Δ ρg = 0.033 m ) (Xu et al., 2006). Specifically, we set 20 rows of the vertical grid blocks near the top boundary to be 0.001 m, ten times smaller than the value estimated by the linear stability analysis. Away from the top boundary, the vertical grid spacing gradually increases to 0.01 m in order to achieve better computational efficiency. Uniform grid spacing of 0.01 m is adopted for both the other two dimensions.

(1)

The equilibrium constant of H2O (KH2O) is obtained from Spycher’s correlation (Spycher et al., 2003) while for all the other species, including CO2 and the relevant impurities, the equilibrium constants are calculated based on temperature, pressure and compositions (Ziabakhsh-Ganji and Kooi, 2012). The Rachford-Rice equation can be solved iteratively by either the Newton’s method or the bisection method. Once the solution β is obtained, the component mole fractions of both phases can be calculated as:

xi =

zi , 1 + β (Ki − 1)

yi = Ki x i

(2) Fig. 1. Mathematical model representation, including the initial and boundary conditions. Figure modified from Li and Jiang (Li and Jiang, 2014).

The obtained dissolved gas compositions (Table 2) are then used as the top boundary conditions. In the base case, the diffusion coefficients 3

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

starting time step (0.01 s) is used to capture the early-time small-scale behavior and to accurately estimate the onset time of convection. In order to achieve better computational efficiency, an automatic time stepping algorithm is also employed (Noels et al., 2002).

Table 1 Simulation parameters. Fluid properties Pressure Temperature Pure water density Density increase from CO2 dissolution Pure water viscosity Salinity Formation properties Permeability Porosity Model parameter Domain length L Domain width W Domain height H

100 bar 45 °C 994.56 kg/m3 10.51 kg/m3 0.5975 × 10−3 Pa·s 0

3. Definition of characteristic parameters Rayleigh number (Ra) is an important non-dimensional parameter to characterize the natural convection in porous media (EmamiMeybodi et al., 2015), which can be expressed as follows,

1.0 × 10−11 m2 0.3

Ra = 1m 1m 1m

100% CO2 90% CO2 + 90% CO2 + 90% CO2 + 90% CO2 + 90% CO2 + 90% CO2 +

10% 10% 10% 10% 10% 10%

Ar CH4 H2S N2 O2 SO2

Dissolved gas composition CO2

Impurity

2.094 × 10−2 1.883 × 10−2 1.871 × 10−2 1.950 × 10−2 1.876 × 10−2 1.880 × 10−2 2.094 × 10−2

– 3.067 × 10−4 2.293 × 10−4 4.642 × 10−3 1.373 × 10−4 2.748 × 10−4 1.018 × 10−2

(7)

In this equation, Δρ is the density difference between gas-saturated formation liquid and gas-free formation liquid at the initial condition, considering the effects of both CO2 and impurity. Rayleigh number compares the strength of the convective mechanism over the diffusion one. The occurrence of density-driven natural convection is expected when the diffusive boundary layer is thick enough and Rayleigh number exceeds a critical value (Emami-Meybodi et al., 2015). The onset time of convection is one of the essential operation issues in assessing the potential of the storage sites. Before the occurrence of density-driven convection, transport of CO2 is controlled by molecular diffusion. The dissolution rate in this period, equal to diffusive flux, declines with time as t−1/2 (Pruess and Zhang, 2008). When convective mixing occurs, the dissolution rate of CO2 increases, which will shorten the dissolving time of mobile CO2 streams into the formation aquifers. Linear stability analyses suggest that the dimensional onset time of convection is related to the properties of the sequestration aquifer as well as the fluid, and can be expressed by the following formula (EnnisKing and Paterson, 2003, 2005; Hassanzadeh et al., 2007; Xu et al., 2006),

Table 2 Injected gas compositions and the corresponding dissolved gas compositions at 45 °C and 100 bar. The units are mole fractions. Injected gas compositions

ΔρgkH ϕμD

The numerical scheme used in the study is based on a fully implicit formulation with Newton-Raphson iteration (Li and Jiang, 2014). Specifically, the governing equations are discretized in space using the integral finite difference method while time is discretized fully implicitly using a first-order backward finite difference scheme. A small

tc = c 0

(ϕμ)2D 1 H2 = c0 2 (Δρgk )2 Ra D

(8) Fig. 2. Time evolution of (a) CO2 dissolution flux, (b) Sherwood number, (c) normalized impurity dissolution flux using CO2 dissolution flux in the 100% CO2 case as the reference, and (d) total dissolved gas inventory for different impurity cases. The grey dotted line in (b) represents the Sherwood number corresponding to pure diffusion with value of 1.0.

4

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

that of pure CO2, making the 10% SO2 case being the only one that experiences a much earlier onset time and higher dissolution rate. Since the velocity in the vertical direction is proportional to density increase (Hassanzadeh et al., 2007), the leading fingers in the SO2 case would be the first to reach the bottom of the fixed model height. Specifically, because the density increase of the 10% SO2 case is two times of that of the pure CO2 case, the leading finger migrates two times faster, making the decay time only half of that in the pure CO2 case (Fig. 2b and Table 3). The maximum Sherwood number in the SO2 case is also distinctively larger, corresponding to stronger convective mixing. The curves of CO2 dissolution flux and Sherwood number of the four cases with the same content of non-condensable impurities, including Ar, CH4, N2, and O2, almost coincide with each other, suggesting that the effects of these four non−CO2 species on the onset and development of convective mixing are comparable. For the 10% H2S case, even though dissolved CO2 concentration at the top boundary is higher than those non-condensable impurity cases (Table 2), the characteristics of the convective mixing process in the H2S case are generally comparable with these non-condensable impurity cases, including density increase, Rayleigh number, onset time of convection, Sherwood number, and decay time of convection (Fig. 2 and Table 3). Although the co-injected SO2 is only 10%, SO2 dissolution rate is generally comparable with CO2 dissolution rate in the pure CO2 case because of its preferential solubility and positive contribution to the density-driven convective mixing. As a result, the total amount of dissolved gases in the 10% SO2 case is significantly enhanced compared with the pure CO2 case (Fig. 2d). On the contrary, the dissolution rate of all the non-condensable impurities is almost two-orders lower because of their lower solubility and negative effects on the density-driven convection. The dissolved quantity of low-solubility impurity is thus negligible and the total dissolved gas inventory in these cases is thus reduced and comparable with each other. Even though both the solubility (Table 2) and dissolution rate (Fig. 2c) of H2S are more than oneorder higher than these non-condensable impurities, the curve of H2S dissolution rate shares similar variation trend with these low-solubility impurities, an implication of similar onset time and decay time. Still, total dissolved gas inventory in the H2S case exceeds that in the pure CO2 case because of its preferential solubility.

There are large discrepancies in the estimation of onset time in the literature according to this equation, which can be attributed to the different methods, initial conditions (the nature and the amplitude of perturbations triggering the convection), and criteria used to define the onset time (Emami-Meybodi et al., 2015; Hassanzadeh et al., 2007). Specially, the numerical constant c0 is case-sensitive, which lies in a rather wide range, from less than 100 to over 5000 (Emami-Meybodi et al., 2015; Pau et al., 2010; Rasmusson et al., 2017). In order to obtain better agreements with experiments, it is suggested that the onset time should be defined as the time at which the dissolution flux begins to grow and/or fingers becomes visible (Thomas et al., 2018). In this study, the critical onset time of convection (tc) is defined as the time when the dissolution rate of CO2 across the top boundary begins to grow, i.e., exceeding the pure diffusive flux. With H2/D as the scale, this dimensional onset time can be further scaled to dimensionless onset time τc. Sherwood number (Sh) is a non-dimensional measure of the relative contribution of convective to diffusive flux during the convective mixing process, which is defined as follows (Emami-Meybodi et al., 2015),

Sh =

Total mixing flux Convective flux + diffusive flux = Diffusive flux Diffusive flux

(9)

It can be used to evaluate and quantify the strength and magnitude of convection in dissolution trapping from the beginning to the end of the dissolution process. For example, it remains unity before the onset of convection. When the convective mixing begins, however, a convective flux starts, and thus Sherwood number starts to exceed unity. During the steady convective mixing period, Sherwood number increases with time since the convective mixing is becoming the primary mass transfer mechanism. After the convective fingers reach the bottom boundary, however, the convective mixing is weakened, which signifies the decay of convection. Both dissolution rate and Sherwood number would decrease accordingly. We define the decay time of convection (denoted as td) as the time when Sherwood number reaches its maximum, i.e., Shmax. This dimensional decay time is also scaled by H2/D to obtain a dimensionless quantity (τd). It should be noted that although this study considers non−CO2 species, CO2 remains the dominant component and the migration of CO2 remains as the focus in the present investigation. Unless otherwise noted, the properties and parameters refer to those of CO2, including the dissolution rate, diffusion coefficient, and Sherwood number.

4.2. Using representative impurity for multiple impurities The properties of the impure streams are related to their compositions, which may affect the occurrence and development of the convective mixing process. Measurement of the composition-dependent parameters are subject to experimental determination and sometimes in need of special materials and equipment. In addition, modeling and simulation of these properties for various mixtures as well as the phase equilibrium and transition could be laborious and time-consuming. It is suggested that the physical properties including phase envelope and density calculation of CO2 mixtures containing air-derived gas impurities, including N2, O2 and Ar, could be represented by that of a single impurity component, such as Ar or N2, at similar concentration level (Wang et al., 2011a). If the effects of CO2 mixtures containing these air-derived impurities on the convective mixing process could be

4. Results and discussions 4.1. Effects of co-injected impurity Fig. 2 shows the time evolution of CO2 dissolution flux, Sherwood number, normalized impurity dissolution flux using CO2 dissolution flux in the 100% CO2 case as the reference, and total dissolved gas inventory for different impurity cases. The corresponding density increase, Rayleigh number, onset time, maximum Sherwood number, and decay time are listed in Table 3. Among all the common impurities, only SO2 would result in larger density increase of the aqueous phase than Table 3 Results for different impurity cases in Table 2. Case 100% CO2 10% Ar 10% CH4 10% H2S 10% N2 10% O2 10% SO2

Δρ (kg/m3) 10.51 9.65 9.13 9.44 9.64 9.46 23.84

Ra 2876 2641 2499 2583 2638 2589 6522

τc

tc (s)

Shmax −4

4

8.17 × 10 8.77 × 104 1.00 × 105 8.92 × 104 9.54 × 104 9.30 × 104 1.51 × 104

1.63 × 10 1.75 × 10−4 2.00 × 10−4 1.78 × 10−4 1.91 × 10−4 1.86 × 10−4 3.02 × 10−5

5

6.71 6.40 6.27 6.46 6.36 6.63 8.73

τd

td (s) 6

2.72 × 10 2.92 × 106 3.34 × 106 3.04 × 106 3.06 × 106 3.05 × 106 1.35 × 106

5.45 × 10−3 5.84 × 10−3 6.68 × 10−3 6.07 × 10−3 6.12 × 10−3 6.10 × 10−3 2.70 × 10−3

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

evaluated using only one impurity component, the modeling and simulation, and laboratory experiments can be substantially simplified. It is for these reasons that we compare the results of single-impurity case with that of multi-impurity (N2, O2, and Ar) mixtures at equivalent impurity molar fractions. N2 is considered as the representative impurity because it is usually the most common impurity in the co-injected CO2 streams and the density increase of the aqueous phase caused by the CO2-N2 mixture is in between the values for the Ar and O2 cases (Table 3). Fig. 3 compares the dissolution rate of CO2 and Sherwood number involving two and three impurities with that of single N2 impurity case, with the corresponding results of the pure CO2 case as references. The curves of these impurity cases overlap pretty well with each other, indicating that the 10% N2 could be used as the representative of the multiple air-derived impurities at the same fraction. The proposed approach of using single N2 impurity to represent multiple N2-O2-Ar impurities seems reasonable, and could be very useful for faster estimation of the effects of these impurities especially from oxy-fuel combustion on the convective mixing process in impure CO2 geological storage. In addition, the reactive and potentially hazardous O2 could also be avoided. Nevertheless, experimental validation should be conducted before applying this approach to more general scenarios.

Table 4 Scenarios for the examination of the effects of gas diffusion coefficients in pure and impure CO2 cases. Case

100% CO2

Impure CO2

Scenario

DCO2 (m2/s)

Dimp (m2/s)

0.5D D 1.5D 2D 3D D-D (Base) 2D-D D-2D 2D-2D

1 × 10−9 2 × 10−9 3 × 10−9 4 × 10−9 6 × 10−9 2 × 10−9 4 × 10−9 2 × 10−9 4 × 10−9

– – – – – 2 × 10−9 2 × 10−9 4 × 10−9 4 × 10−9

a result, the time needed for the occurrence of convection increases with increasing CO2 diffusivity. Fig. 5 illustrates the relations between different characteristic parameters for these scenarios in pure CO2 case. Fig. 5a demonstrates that there is a linear relation between the dimensional onset time and CO2 diffusion coefficient, which agrees with both previous theoretical linear analyses (Eq. (8)) and numerical simulation results (Pau et al., 2010). Fig. 5b shows that the dimensionless onset time scales with Ra as (Ra)−2.020, which agrees pretty well with τc ∼ (Ra)-2 in the previous investigations (Emami-Meybodi et al., 2015). This agreement suggests that the definition of the onset time in the present study is reasonable and comparable with theoretical analyses. As shown in Fig. 5c, the maximum Sherwood number may be estimated based on the Rayleigh number by 0.326Ra0.378. In previous investigations, however, it was suggested that Shmax = 0.105Ra0.5 (Hassanzadeh et al., 2007) or Shmax = 0.115Ra0.77 (Moghaddam et al., 2015). Fig. 5d shows that the dimensionless decay time is inversely proportional to Rayleigh number as Ra-0.922. Similarly, previous investigations implied that dimensionless time at which the maximum Sherwood number occurs was 100Ra6/5 (Hassanzadeh et al., 2007) or 0.91Ra-0.45 (Moghaddam et al., 2015). The differences in the pre-factors and/or exponents between our expressions and previous studies may result from the different definitions of Sherwood number. Although most studies define Sherwood number based on the ratio of total mixing to mixing achieved by diffusion, different investigations deal with the definition in slightly different ways. For example, Sherwood number stands for the ratio of the sum of convective and diffusive flux to the diffusive flux across the interface in the present study while some researchers obtained it using mathematical expression derived from their specific experiments (EmamiMeybodi et al., 2015; Moghaddam et al., 2015). Nevertheless, the general variation trend and pattern of these scaling relationships are similar. For impure CO2 case, three different pairs of diffusion coefficients

4.3. Effects of diffusion coefficients The density increase of the formation fluid caused by gas dissolution is the driving force of the convective mixing process. Meanwhile, the strength and characteristics of the convective mixing are related to Rayleigh number which is dependent on diffusion coefficient in addition to the density increase (Eq. (7)). However, available diffusion coefficients of gases in water measured from different methods are quite different from each other (Winkelmann, 2007). Direct measurements of gas diffusivity in the formation fluids are very complex and require special equipment (Mahmoodpour et al., 2018). In this section, we intend to investigate the effects of diffusion coefficients of different gases on the characteristics of the convective mixing process. The pure CO2 case is first examined and CO2 diffusion coefficient is varied from 1 × 10−9 to 6 × 10−9 m2/s (Table 4). Fig. 4 illustrates the variation of CO2 dissolution flux and Sherwood number for these different scenarios in pure CO2 case. Rayleigh number, onset time, maximum Sherwood number, and decay time of convection in response to different CO2 diffusivity are shown in Table 5. Before the occurrence of convection, CO2 dissolution rate in the early diffusive period, i.e., pure diffusive flux, increases with increasing CO2 diffusivity as predicted (Emami-Meybodi et al., 2015). Sherwood number (Fig. 4b) and Rayleigh number (Table 5), however, decreases with increasing CO2 diffusion coefficient, indicating a weakened convective mixing process. As

Fig. 3. Time evolution of (a) CO2 dissolution flux and (b) Sherwood number for different cases. The grey dotted line in (b) represents the Sherwood number corresponding to pure diffusion with value of 1.0. 6

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

Fig. 4. Variation of (a) CO2 dissolution flux and (b) Sherwood number with different CO2 diffusion coefficients in pure CO2 case. The grey dotted line in (b) represents the Sherwood number corresponding to pure diffusion with value of 1.0.

are examined for each impurity and compared with the base case in Section 4.1 where both DCO2 and Dimp are 2 × 10−9 m2/s (Table 4). Fig. 6 shows the time evolution of CO2 dissolution flux, impurity flux as well as Sherwood number for the 10% Ar and 10% SO2 case. The variation trends of the other impurity cases are similar as that of the 10% Ar case and thus are not presented. The resulted characteristic parameters are listed in Table 5. For the CO2-imurity cases, the general variation trend of CO2 dissolution rate resulted from doubling the diffusion coefficient of CO2 alone (D-D vs 2D-D and D-2D vs 2D-2D) is similar as that in pure CO2 case (D-100% CO2 vs 2D-100% CO2). However, increasing CO2 diffusivity alone decreases impurity dissolution rate and Sherwood number on the whole for all the investigated scenarios. Specifically, the maximum Sherwood number decreases significantly by doubling CO2

diffusivity alone. For instance, Shmax in the 10% Ar case decreases by more than 20% when CO2 diffusivity is doubled. Except for the 10% SO2 case, dimensional onset time of convection is more than two times larger as a result of doubled CO2 diffusivity alone while the dimensionless figure is over fourfold. In the 10% SO2 case, the increase of the onset time is relatively smaller than other cases. Except for the 10% SO2 case, doubling the diffusion coefficient of impurity alone (D-D vs D-2D or 2D-D vs 2D-2D) increases impurity dissolution rate but does not have distinct effects on CO2 dissolution rate or Sherwood number. In fact, varying diffusion coefficient of these impurities is not likely to have influential effects on the characteristics of convective mixing process, especially the onset time and the maximum Sherwood number Shmax (Table 5). For the 10% SO2 case, however, with increasing diffusion coefficient of SO2 impurity alone,

Table 5 Rayleigh number, onset time, maximum Sherwood number and decay time of convection for different diffusivity scenarios in pure and impure CO2 cases. As defined in Table 4, for 100% CO2 cases, D refers to the base case with CO2 diffusivity at 2 × 10−9 m2/s while the numbers in front of D in other scenarios stand for its multiplication factor. For impure CO2 cases, D-D means that the diffusvities of CO2 and impurity are equal as 2 × 10−9 m2/s, with the first D for CO2 and the second for impurity. In this sense, 2D-D means doubling the CO2 diffusivity only, D-2D means doubling the impurity diffusivity only, and 2D-2D means doubling both CO2 and impurity diffusivities to 4 × 10−9 m2/s. Case

100% CO2

10% Ar

10% CH4

10% H2S

10% N2

10% O2

10% SO2

Scenario 0.5D D 1.5D 2D 3D D-D D-2D 2D-D 2D-2D D-D D-2D 2D-D 2D-2D D-D D-2D 2D-D 2D-2D D-D D-2D 2D-D 2D-2D D-D D-2D 2D-D 2D-2D D-D D-2D 2D-D 2D-2D

Ra 5753 2876 1917 1438 959 2641 2641 1320 1320 2499 2499 1249 1249 2583 2583 1292 1292 2638 2638 1319 1319 2589 2589 1295 1295 6522 6522 3261 3261

τc

tc (s)

Shmax −5

4

3.94 × 10 8.17 × 104 1.11 × 105 1.60 × 105 2.52 × 105 8.77 × 104 9.15 × 104 2.00 × 105 2.00 × 105 1.00 × 105 1.00 × 105 2.20 × 105 2.19 × 105 8.92 × 104 8.92 × 104 2.04 × 105 2.00 × 105 9.54 × 104 8.85 × 104 1.93 × 105 2.00 × 105 9.30 × 104 9.75 × 104 2.07 × 105 2.00 × 105 1.51 × 104 2.04 × 104 2.10 × 104 2.93 × 104

3.94 × 10 1.63 × 10−4 3.33 × 10−4 6.40 × 10−4 1.51 × 10−3 1.75 × 10−4 1.83 × 10−4 8.00 × 10−4 8.00 × 10−4 2.00 × 10−4 2.00 × 10−4 8.82 × 10−4 8.77 × 10−4 1.78 × 10−4 1.78 × 10−4 8.18 × 10−4 8.00 × 10−4 1.91 × 10−4 1.77 × 10−4 7.71 × 10−4 8.00 × 10−4 1.86 × 10−4 1.95 × 10−4 8.27 × 10−4 8.00 × 10−4 3.02 × 10−5 4.09 × 10−5 8.42 × 10−5 1.17 × 10−4

7

8.52 6.71 5.84 4.90 4.33 6.40 6.48 4.82 4.81 6.27 6.31 4.94 4.92 6.46 6.53 5.01 4.64 6.36 6.57 4.76 4.82 6.63 6.58 4.74 4.77 8.73 8.70 7.08 6.86

τd

td (s) 6

3.08 × 10 2.72 × 106 2.76 × 106 2.28 × 106 2.88 × 106 2.92 × 106 3.20 × 106 2.72 × 106 3.18 × 106 3.34 × 106 3.29 × 106 2.92 × 106 3.04 × 106 3.04 × 106 3.21 × 106 2.92 × 106 3.09 × 106 3.06 × 106 3.08 × 106 2.83 × 106 3.16 × 106 3.05 × 106 3.12 × 106 3.20 × 106 3.16 × 106 1.35 × 106 1.22 × 106 1.26 × 106 1.21 × 106

3.08 × 10−3 5.45 × 10−3 8.28 × 10−3 9.11 × 10−3 1.73 × 10−2 5.84 × 10−3 6.40 × 10−3 1.09 × 10−2 1.27 × 10−2 6.68 × 10−3 6.58 × 10−3 1.17 × 10−2 1.22 × 10−2 6.07 × 10−3 6.43 × 10−3 1.17 × 10−2 1.24 × 10−2 6.12 × 10−3 6.15 × 10−3 1.13 × 10−2 1.26 × 10−2 6.10 × 10−3 6.25 × 10−3 1.28 × 10−2 1.26 × 10−2 2.70 × 10−3 2.43 × 10−3 5.03 × 10−3 4.85 × 10−3

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

Fig. 5. Relations between different characteristic parameters in pure CO2 case, (a) dimensional onset time and CO2 diffusion coefficient, (b) dimensionless onset time and Rayleigh number, (c) maximum Sherwood number and Rayleigh number, (d) dimensionless decay time and Rayleigh number.

with SO2, SO2 diffusivity should be taken into consideration in the calculation of effective Rayleigh number to obtain more reasonable predictions, especially for scenarios with unequal CO2 and SO2 diffusivity. The results of 10% SO2 case with non-equal CO2 and SO2 diffusivity are also excluded to obtain the relationships between maximum Sherwood number and Rayleigh number as well as dimensionless decay time and Rayleigh number, as shown in Fig. 7c and d, respectively,

SO2 dissolution rate increases while CO2 dissolution rate decreases on a whole. Specifically, onset time tc increases significantly with increasing SO2 diffusivity, indicating weakened convective mixing of CO2-rich aqueous phase, even though both the Rayleigh number as well as the Sherwood number do not exhibit distinguishable change. This is because these two characteristic numbers are calculated based on the properties of CO2, not accounting for the effects of the non−CO2 species (Section 3). Fig. 7 demonstrates the relationships between different characteristic parameters when taking into account both the effects of density increase and diffusion coefficients of CO2 and impurity (Table 5). Compared with results in Fig. 5 for pure CO2 case only, it is suggested that the existence of impurity alters the scaling relationship of these characteristic parameters to different extents. Dimensionless onset time in both pure CO2 and impure CO2 cases can be scaled with Rayleigh number as τc∼Ra−2.031, as shown in Fig. 7a. Closer examination shows that the largest deviations come from the two scenarios of non-equal CO2 and SO2 diffusivity (D-2D and 2D-D) in 10% SO2 case. As mentioned above, the Rayleigh number calculated in the present study is related to the apparent CO2 diffusivity only, neglecting the diffusivity of co-injected non−CO2 species. However, practical CO2 diffusivity was suggested to be affected by impurity diffusivity (Mahmoodpour et al., 2018). It can then be deduced that the 10% SO2 may affect practical CO2 diffusivity and the effective Rayleigh number of the system while the other impurities do not have significant effects. The results of these two scenarios in the 10% SO2 case are then removed, and the new fitting equation is shown in Fig. 7b as:

τc = 2147.83Ra −2.065

Shmax = 0.328Ra0.377

(11)

τd = 9.84Ra −0.937

(12)

Different impurities would result in different density increases of the formation liquid, and the case would be more complicated taking into account different diffusivities of CO2 and non−CO2 species. The scaling relationships showing in Fig. 7 imply that the onset, strength, and decay of the convective mixing process could be estimated by the Rayleigh number calculated by the dominating CO2 properties and formation properties for most CO2-impurity scenarios. However, special attention should be paid to the cases involving SO2. Apart from the features of high-solubility and high-density-increase of the aqueous phase, the effects of SO2 on the practical CO2 diffusivity and thus the effective Rayleigh number may not be neglected as other common impurities. 4.4. Discussion Our model and numerical simulations are based on several assumptions. To begin with, the porous medium is simplified to an isotropic and homogenous one with a small random permeability modification ( ± 4%) which triggers the onset of convection. However, real geological formations are far more complicated and are usually anisotropic and heterogeneous. It is suggested that the permeability heterogeneity alone may have a remarkable influence on the development and growth of the convective fingers (Emami-Meybodi et al., 2015). For example, the heterogeneous formations may experience much earlier

(10) −2

The scaling relation is similar as τc∼Ra in pure CO2 case found in a previous study (Emami-Meybodi et al., 2015). Compared with the fitting relations in Fig. 5b considering only CO2 diffusivity in pure CO2 case, it is implied that the diffusivity of non−CO2 species (excluding SO2) has an insignificant influence in the estimation of dimensionless onset time based on Rayleigh number. However, for cases co-injected 8

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

Fig. 6. Variation of CO2 dissolution rate, impurity dissolution rate, and Sherwood number for 10% Ar case (left column) and 10% SO2 case (right column). Other impurity cases share similar trends with that of 10% Ar case and are not presented here. As defined in Table 4, D100% CO2 and 2D-100% CO2 refer to the pure CO2 scenarios with CO2 diffusivities at 2 × 10−9 m2/s and 4 × 10−9 m2/s, respectively. For impure CO2 scenarios, D-D means that the diffusvities of CO2 and impurity are equal as 2 × 10−9 m2/s, with the first D for CO2 and the second for impurity. In this sense, 2D-D means doubling the CO2 diffusivity only, D-2D means doubling the impurity diffusivity only, and 2D-2D means doubling both CO2 and impurity diffusivities to 4 × 10−9 m2/s.

the capillary transition zone alone has much stronger effects on the onset time and mass transfer rate than the co-injected impurities. In addition, the effects of geochemical reactions, aquifer dipping, temperature gradients, and salinity have also been neglected. For smallscale problems as in the present study, isothermal conditions should not cause any issues. However, for practical field-scale projects, temperature gradients in the reservoirs would definitely affect the solubility and dissolution of CO2 as well as non−CO2 species. Besides, the salinity in the deep saline aquifers reduces gas solubility in the aqueous phase and thus the density change, which is not favorable for the development of the density-driven fingering process and the comparison of different impurities. As a result, we have neglected the salinity in this investigation on the effects of impurities. However, the level of salinity was implied to decrease the effective CO2 diffusion coefficient (Mahmoodpour et al., 2018). In addition, the impure cases may be more sensitive to salinity than pure CO2 cases. It is still not clear what the coupled effects of salinity and impurities on the effective CO2 diffusion coefficient will be. It thus seems necessary to accurately measure or calculate the effective CO2 diffusion coefficient in brines with different salinity and brine compositions for practical impure CO2 sequestration investigations (Mahmoodpour et al., 2019). These limitations should be

onset of convection than the equivalent homogenous formations with same effective properties (e.g., Green and Ennis-King, 2010a, 2010b, 2014, 2018). Besides, heterogeneity was implied to affect the growth of the convective fingers and future upscaling of fine-scale models to fieldscale models. Above all, compared with the real-world reservoirs, the permeability heterogeneity of most existing studies is small or moderate, while larger heterogeneity may have stronger effects on the convective mixing process. Since this study is mainly focused on the fluid flow and transfer in the aqueous phase, the more complicated two-phase flow that may involve potential phase transition as well as capillary transition zone and relative permeability is approximated by a single-phase flow maintaining gas-saturated top boundary, which might have an influence in the evolution of the convective mixing mechanism. For example, it is suggested that neglecting the capillary transition zone would significantly underestimate the contribution of the convective mixing to the security of CO2 geological sequestration (Elenius et al., 2012). Specifically, when taking into consideration the cross-flow between the capillary transition zone and the diffusive boundary layer, the onset time of instability might be reduced by a factor of five while the dissolution rate could be enhanced up to four times. It seems that 9

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

Fig. 7. Relations between different characteristic parameters for both pure CO2 and impure CO2 cases, (a) dimensionless onset time and Rayleigh number for all the results in Table 5, (b) dimensionless onset time and Rayleigh number excluding D-2D and 2D-D scenarios in 10% SO2 case, (c) maximum Sherwood number and Rayleigh number excluding D-2D and 2DD scenarios in 10% SO2 case, (d) dimensionless decay time and Rayleigh number excluding D2D and 2D-D scenarios in 10% SO2 case.

of density-driven convection. More importantly, scaling relationships are presented to characterize the onset, strength, and decay of the convective mixing process in impure CO2 geological storage. These relations can be used in selection of the best potential aquifers among different candidate reservoirs and the implementation of large-scale storage in deep saline aquifers. Special attention should be paid to cases involving SO2, however, because there would be pretty distinct deviations from the fitting relations, which may be because the Rayleigh number calculated in the present study is based on the apparent CO2 diffusivity only, neglecting the effects of the co-injected impurity. The neglect of the effects of other common impurities investigated seems unlikely to make much difference on CO2 diffusivity and thus the predictions but the existence of SO2 is supposed to influence the practical CO2 diffusivity based on our investigations. When SO2 is involved, it seems necessary to measure or estimate the effective CO2 diffusivity to obtain more reasonable predictions of effective Rayleigh number and scaling relations.

kept in mind when applying the analyses and models to other systems or specific reservoirs. Further studies are needed to take into consideration of anisotropy, heterogeneity, capillary transition zone, temperature gradient, salinity, and multi-phase effects on the scaling behavior of convective mixing in impure geological storage. Nevertheless, our results demonstrated in this investigation could give insight and a better understanding of convective mixing process in impure CO2 geological storage as a simplified case. The scaling relations can provide reference for the selection of appropriate candidate reservoirs and the implementation of large-scale storage in deep saline aquifers.

5. Conclusions The density-driven convective mixing has been demonstrated to significantly improve CO2 dissolution rate into the formation brine and plays an important role in efficient and secure long-term geological storage. Most practical CO2 geological projects involve in several non−CO2 species, such as N2, O2, Ar, H2S, and SO2. This study investigates and compares the effects of different impurities on the convective mixing process. Comparing with the pure CO2 case, SO2 is the only impurity that increases the density increase of the aqueous phase, accelerates the onset of convection, and increases the Rayleigh number and the maximum Sherwood number (Table 3). The dissolution rate as well as the total dissolved gas inventory in the CO2-SO2 case also significantly increase compared with all the other cases (Fig. 2). Meanwhile, all the other common impurities investigated would lead to delayed occurrence of convection and reduced CO2 dissolution rate. Specially, it is suggested that single impurity such as N2 or Ar might be used as representative for multiple air-derived species in estimating the evolution of the convective mixing process. Apart from the effects on density increase, the co-injected species are supposed to affect effective CO2 diffusivity, altering the Rayleigh number and thus the evolution of the convective mixing process. However, our simulation results indicate that except for SO2, it seems unlikely that the diffusion coefficients of all the other impurities examined would have distinct effects on the occurrence and development

CRediT authorship contribution statement Didi Li: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing - original draft. Xi Jiang: Formal analysis, Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This study is supported by the National Natural Science Foundation of China (Grant No. 41807191) and the Natural Science Foundation of Guangdong Province, China (Grant No. 2018A030310524). DL would like to acknowledge the colleagues and friends at the Indiana 10

International Journal of Greenhouse Gas Control 96 (2020) 103015

D. Li and X. Jiang

Geological & Water Survey for providing the perfect working environment. DL also acknowledges financial support from the visiting scholar program of the China Scholarship Council (No. 201808440020).

data and theoretical models. Appl. Energy 88 (11), 3567–3579. Li, D., Jiang, X., Meng, Q., Xie, Q., 2015. Numerical analyses of the effects of nitrogen on the dissolution trapping mechanism of carbon dioxide geological storage. Comput. Fluids 114, 1–11. Liu, Y., Ding, T., Yu, B., Yang, Y., 2017. Convective dissolution analysis of long-term storage of acid gas in saline aquifers. Energy Procedia 114, 3417–3431. Mahmoodpour, S., Rostami, B., Kachouei, S., 2017. Impure CO2 dissolution during geosequestration in saline aquifers with background flow. 79th EAGE Conference and Exhibition 2017. Mahmoodpour, S., Rostami, B., Emami-Meybodi, H., 2018. Onset of convection controlled by N2 impurity during CO2 storage in saline aquifers. Int. J. Greenh. Gas Control. 79, 234–247. Mahmoodpour, S., Rostami, B., Soltanian, M.R., Amooie, M.A., 2019. Effect of brine composition on the onset of convection during CO2 dissolution in brine. Comput. Geosci. 124, 1–13. Mao, S., Duan, Z., 2006. A thermodynamic model for calculating nitrogen solubility, gas phase composition and density of the N2–H2O–NaCl system. Fluid Phase Equilib. 248 (2), 103–114. Metz, B., Davidson, O., De Coninck, H., Loos, M., Meyer, L., 2005. IPCC Special Report on Carbon Dioxide Capture and Storage. Prepared by Working Group III of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. Moghaddam, R.N., Rostami, B., Pourafshary, P., 2015. Scaling analysis of the convective mixing in porous media for geological storage of CO2: an experimental approach. Chem. Eng. Commun. 202 (6), 815–822. Noels, L., Stainier, L., Ponthot, J.P., Bonini, J., 2002. Automatic time stepping algorithms for implicit numerical simulations of non-linear dynamics. Adv. Eng. Softw. 33 (710), 589–603. Pau, G.S.H., Bell, J.B., Pruess, K., Almgren, A.S., Lijewski, M.J., Zhang, K., 2010. Highresolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers. Adv. Water Resour. 33 (4), 443–455. Pruess, K., Zhang, K., 2008. Numerical Modeling Studies of the Dissolution-diffusionconvection Process. Technical Report LBNL-1243E. Lawrence Berkeley National Laboratory, California. Raad, S.M.J., Hassanzadeh, H., 2016. Does impure CO2 impede or accelerate the onset of convective mixing in geological storage? Int. J. Greenh. Gas Control. 54, 250–257. Raad, S.M.J., Hassanzadeh, H., 2017. Prospect for storage of impure carbon dioxide streams in deep saline aquifers—a convective dissolution perspective. Int. J. Greenh. Gas Control. 63, 350–355. Raad, S.M.J., Hassanzadeh, H., 2018. Comments on the paper “effect of impurities on the onset and the growth of gravitational instabilities in a geological CO2 storage process: linear and nonlinear analyses” M.C. Kim, K.H. Song (2017). Chem. Eng. Sci. 192, 613–618. Rachford, H.H., Rice, J.D., 1952. Procedure for use of electrical digital computers in calculating flash vaporization hydrocarbon equilibrium. J. Pet. Technol. 4, 19–20. Rasmusson, M., Fagerlund, F., Rasmusson, K., Tsang, Y., Niemi, A., 2017. Refractive-lighttransmission technique applied to density-driven convective mixing in porous media with implications for geological CO2 storage. Water Resour. Res. 53 (11), 8760–8780. Spycher, N., Pruess, K., Ennis-King, J., 2003. CO2-H2O mixtures in the geological sequestration of CO2. I. Assessment and calculation of mutual solubilities from 12 to 100 °C and up to 600 bar. Geochim. Cosmochim. Acta 67 (16), 3015–3031. Taheri, A., Torsaeter, O., Lindeberg, E., Hadia, N.J., Wessel-Berg, D., 2018. Qualitative and quantitative experimental study of convective mixing process during storage of CO2 in heterogeneous saline aquifers. Int. J. Greenh. Gas Control. 71, 212–226. Talman, S., 2015. Subsurface geochemical fate and effects of impurities contained in a CO2 stream injected into a deep saline aquifer: what is known. Int. J. Greenh. Gas Control. 40, 267–291. Thomas, C., Dehaeck, S., De Wit, A., 2018. Convective dissolution of CO2 in water and salt solutions. Int. J. Greenh. Gas Control. 72, 105–116. Wang, J., Ryan, D., Anthony, E.J., Wigston, A., 2011a. Effects of Impurities on Geological Storage of CO2. IEAGHG. Wang, J., Ryan, D., Anthony, E.J., Wildgust, N., Aiken, T., 2011b. Effects of impurities on CO2 transport, injection and storage. 10th International Conference on Greenhouse Gas Control Technologies 4, 3071–3078. Winkelmann, J., 2007. Diffusion in Gases, Liquids and Electrolytes: Gases in Gases, Liquids and Their Mixtures. Springer-Verlag, Berlin Heidelberg, pp. 15. Xie, Y., Simmons, C.T., Werner, A.D., 2011. Speed of free convective fingering in porous media. Water Resour. Res. 47 (11), W11501. Xu, X., Chen, S., Zhang, D., 2006. Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv. Water Resour. 29 (3), 397–407. Zhang, W., Li, Y., Omambia, A.N., 2011. Reactive transport modeling of effects of convective mixing on long-term CO2 geological storage in deep saline formations. Int. J. Greenh. Gas Control. 5 (2), 241–256. Ziabakhsh-Ganji, Z., Kooi, H., 2012. An Equation of State for thermodynamic equilibrium of gas mixtures and brines to allow simulation of the effects of impurities in subsurface CO2 storage. Int. J. Greenh. Gas Control. 11, S21–S34.

Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.ijggc.2020.103015. References Alcalde, J., Flude, S., Wilkinson, M., Johnson, G., Edlmann, K., Bond, C.E., Scott, V., Gilfillan, S.M.V., Ogaya, X., Haszeldine, R.S., 2018. Estimating geological CO2 storage security to deliver on climate mitigation. Nat. Commun. 9 (1), 2201. Bachu, S., 2008. CO2 storage in geological media: role, means, status and barriers to deployment. Prog. Energy Combust. Sci. 34 (2), 254–273. Chadwick, A., Arts, R., Bernstone, C., May, F., Thibeau, S., Zweigel, P., 2008. Best practice for the storage of CO2 in saline aquifers–observations and guidelines from the SACS and CO2STORE projects. NERC. Crandell, L.E., Ellis, B.R., Peters, C.A., 2010. Dissolution potential of SO2 co-injected with CO2 in geologic sequestration. Environ. Sci. Technol. 44 (1), 349–355. Elenius, M.T., Nordbotten, J.M., Kalisch, H., 2012. Effects of a capillary transition zone on the stability of a diffusive boundary layer. IMA J. Appl. Math. 77 (6), 771–787. Emami-Meybodi, H., Hassanzadeh, H., Green, C.P., Ennis-King, J., 2015. Convective dissolution of CO2 in saline aquifers: progress in modeling and experiments. Int. J. Greenh. Gas Control. 40, 238–266. Ennis-King, J., Paterson, L., 2003. In: Gale, J., Kaya, Y. (Eds.), Rate of Dissolution Due to Convective Mixing in the Underground Storage of Carbon Dioxide, pp. 1653–1656. Ennis-King, J., Paterson, L., 2005. Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations. Spe J. 10 (3), 349–356. Ennis-King, J., Preston, I., Paterson, L., 2005. Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids 17 (8), 084107. Farajzadeh, R., Salimi, H., Zitha, P.L.J., Bruining, H., 2007. Numerical simulation of density-driven natural convection in porous media with application for CO2 injection projects. Int. J. Heat Mass Transf. 50, 5054–5064. Gale, J., 2009. Impure thoughts. Int. J. Greenh. Gas Control. 3 (1), 1–2. Garcia, J.E., 2001. Density of Aqueous Solutions of CO2 Lawrence Berkeley National Laboratory Technical Report LBNL-49023. Green, C.P., Ennis-King, J., 2010a. Effect of vertical heterogeneity on long-term migration of CO2 in saline formations. Transp. Porous Media 82 (1), 31–47. Green, C.P., Ennis-King, J., 2010b. Vertical permeability distribution of reservoirs with impermeable barriers. Transp. Porous Media 83 (3), 525–539. Green, C.P., Ennis-King, J., 2014. Steady dissolution rate due to convective mixing in anisotropic porous media. Adv. Water Resour. 73, 65–73. Green, C.P., Ennis-King, J., 2018. Steady flux regime during convective mixing in threedimensional heterogeneous porous media. Fluids 3, 58. Hassanzadeh, H., Pooladi-Darvish, M., Keith, D.W., 2007. Scaling behavior of convective mixing, with application to geological storage of CO2. Aiche J. 53 (5), 1121–1131. Jacquemet, N., Pironon, J., Saint-Marc, J., 2007. Mineralogical changes of a well cement in various H2S-CO2(-brine) fluids at high pressure and temperature. Environ. Sci. Technol. 2008 (42), 282–288. Ji, X., Zhu, C., 2013. Predicting possible effects of H2S impurity on CO2 transportation and geological storage. Environ. Sci. Technol. 47 (1), 55–62. Jiang, X., 2011. A review of physical modelling and numerical simulation of long-term geological storage of CO2. Appl. Energy 88 (11), 3557–3566. Johnson, J.W., Oelkers, E.H., Helgeson, H.C., 1992. SUPCRT92: a software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000°C. Comput. Geosci. 18 (7), 899–947. Kim, M.C., Song, K.H., 2017. Effect of impurities on the onset and growth of gravitational instabilities in a geological CO2 storage process: linear and nonlinear analyses. Chem. Eng. Sci. 174, 426–444. Kim, M.C., Song, K.H., 2019. Responses to the comment on the paper “Effect of impurities on the onset and the growth of gravitational instabilities in a geological CO2 storage process: linear and nonlinear analyses” by M.C. Kim and K.H. Song. Chem. Eng. Sci. 193, 184–187. Leibovici, C.F., Neoschil, J., 1995. A solution of Rachford-Rice equations for multiphase systems. Fluid Phase Equilib. 112, 217–221. Li, D., Jiang, X., 2014. A numerical study of the impurity effects of nitrogen and sulfur dioxide on the solubility trapping of carbon dioxide geological storage. Appl. Energy 128, 60–74. Li, H., Jakobsen, J.P., Wilhelmsen, Ø., Yan, J., 2011. PVTxy properties of CO2 mixtures relevant for CO2 capture, transport and storage: review of available experimental

11