- Email: [email protected]
Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions
Accepted Manuscript
Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions Maria Rasmusson , Fritjof Fagerlund , Yvonne Tsang , Kristina Rasmusson , Auli Niemi PII: DOI: Reference:
S0309-1708(15)00195-5 10.1016/j.advwatres.2015.08.009 ADWR 2446
To appear in:
Advances in Water Resources
Received date: Revised date: Accepted date:
27 October 2014 8 August 2015 11 August 2015
Please cite this article as: Maria Rasmusson , Fritjof Fagerlund , Yvonne Tsang , Kristina Rasmusson , Auli Niemi , Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions, Advances in Water Resources (2015), doi: 10.1016/j.advwatres.2015.08.009
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights We develop an instability onset, onset time and steady convective flux-analysis tool. A systematic analysis of prerequisites for enhanced solubility trapping is presented. Indirect thermal effects on enhanced solubility trapping are non-negligible. Shallow storage depths favor density-driven instability onset and short onset times. Salinity determines the storage depth at which the largest convective fluxes occur.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Title: Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions
CR IP T
Authors: Maria Rasmussona, Fritjof Fagerlunda, Yvonne Tsangb, Kristina Rasmussona, Auli
Corresponding author: M. Rasmussona
a
M
Affiliations:
AN US
Niemia
Uppsala University – Department of Earth Sciences, Villavagen 16, SE-75236 Uppsala,
ED
Sweden.
PT
+46184717162, [email protected]
CE
+46184717166, [email protected]
AC
+46184717162, [email protected] +46184712263, [email protected] fax number +4618551124 (Department of Earth Sciences) b
Lawrence Berkeley National Laboratory, University of California, Earth Sciences Division, 1
Cyclotron Road, Berkeley, CA 94720-8126, USA. +15104867047, [email protected]
ACCEPTED MANUSCRIPT Abstract Direct atmospheric greenhouse gas emissions can be greatly reduced by CO2 sequestration in deep saline aquifers. One of the most secure and important mechanisms of CO2 trapping over large time scales is solubility trapping. In addition, the CO2 dissolution rate is greatly enhanced if density-driven convective mixing occurs. We present a systematic analysis of the
CR IP T
prerequisites for density-driven instability and convective mixing over the broad temperature, pressure, salinity and permeability conditions that are found in geological CO2 storage. The onset of instability (Rayleigh-Darcy number, 𝑅𝑎), the onset time of instability and the steady
AN US
convective flux are comprehensively calculated using a newly developed analysis tool that accounts for the thermodynamic and salinity dependence on solutally and thermally induced density change, viscosity, molecular and thermal diffusivity. Additionally, the relative influences of field characteristics are analysed through local and global sensitivity analyses.
M
The results help to elucidate the trends of the 𝑅𝑎, onset time of instability and steady
ED
convective flux under field conditions. The impacts of storage depth and basin type (geothermal gradient) are also explored and the conditions that favour or hinder enhanced
PT
solubility trapping are identified. Contrary to previous studies, we conclude that the geothermal gradient has a non-negligible effect on density-driven instability and convective
CE
mixing when considering both direct and indirect thermal effects because cold basin
AC
conditions, for instance, render higher 𝑅𝑎 compared to warm basin conditions. We also show that the largest 𝑅𝑎 is obtained for conditions that correspond to relatively shallow depths, measuring approximately 800 m, indicating that CO2 storage at such depths favours the onset of density-driven instability and reduces onset times. However, shallow depths do not necessarily provide conditions that generate the largest steady convective fluxes; the salinity determines the storage depth at which the largest steady convective fluxes occur.
ACCEPTED MANUSCRIPT Furthermore, we present a straight-forward and efficient procedure to estimate site-specific solutal 𝑅𝑎 that accounts for thermodynamic and salinity dependence.
Concentration [𝑘𝑔𝐶𝑂2 ⁄𝑚3 ]
𝑐
Specific heat capacity [ 𝐽 ⁄(𝑘𝑔 ∙ 𝐾)]
𝑐0
Numerically derived constant [-]
𝐷
Molecular diffusivity [𝑚2 ⁄𝑠]
EE
Mean elementary effect
𝐹
Steady convective flux [𝑘𝑔𝐶𝑂2 ⁄(𝑚2 ∙ 𝑠)]
𝐹∗
Steady convective flux [-]
𝑔
Magnitude of gravitational acceleration [𝑚⁄𝑠 2 ]
𝐻
Layer thickness [𝑚]
𝑘
Permeability [𝑚2 ]
M
ED
PT
CE
Molecular weight [𝑔⁄𝑚𝑜𝑙 ]
AC
̃ 𝑀
AN US
𝐶
CR IP T
Nomenclature
𝑁
Buoyancy ratio [-]
P
Pressure [𝑀𝑃𝑎]
r
Number of paths in Morris OAT method
𝑅𝑎
Rayleigh-Darcy number [-]
ACCEPTED MANUSCRIPT Temperature [˚C (if not otherwise specified)]
𝑡𝑐
Onset time of instability [s]
𝑆
S-number [-]
Salinity
Salinity [𝑚𝑜𝑙𝑁𝑎𝐶𝑙 ⁄𝑘𝑔𝑤 ]
SD
Standard deviation
𝑉̅
Molar volume [𝑐𝑚3 ⁄𝑚𝑜𝑙 ]
𝑥
Solubility [𝑚𝑜𝑙𝐶𝑂2 ⁄𝑘𝑔𝑤 ]
AN US
CR IP T
T
Greek letters
Induced density change [𝑘𝑔⁄𝑚3 ] or [%]
𝛻𝑇𝑔𝑒𝑜
Geothermal gradient [°𝐶 ⁄𝑘𝑚]
𝜃
Empirical constant
𝜅
Thermal diffusivity [𝑚2 ⁄𝑠]
ED
PT
CE
Thermal conductivity of medium [𝑊/(𝑚 ∙ 𝐾)]
AC
𝜆
M
∆𝜌
𝜇
Viscosity [𝑘𝑔⁄(𝑚 ∙ 𝑠)]
𝜌
Density [𝑘𝑔⁄𝑚3 ]
𝜑
Porosity [-]
ACCEPTED MANUSCRIPT
Critical
𝑓
Formation fluid
𝑆
Solutal
𝑠𝑎𝑡
Saturation
𝑠𝑢𝑟𝑓
Surface (temperature)
𝑇
Thermal
𝑤
Water
%
Percentage
1. Introduction
ED
M
AN US
𝑐
CR IP T
Subscripts
PT
CO2 storage in sedimentary basins that contain deep saline aquifers reduces the direct release of greenhouse gases into the atmosphere. The method entails injecting CO2 at ~1-2 km depth
CE
into formations of sufficient porosity, permeability and injectivity. Depending on the
AC
thermodynamic conditions, the injected CO2-rich phase is in a gaseous, liquid or supercritical state. The injected CO2 is retained in the formation by several mechanisms: physical (structural, stratigraphic and hydrodynamic), residual, solubility and mineral trapping [1]. Over large time scales, solubility trapping is one of the most important trapping mechanisms, providing more secure storage compared to free phase CO2 containment [1]. However, CO2 solubility in the aqueous phase (brine) depends on the temperature, pressure and salinity, and the same applies to the accompanying density change in the aqueous phase. Density increases
ACCEPTED MANUSCRIPT of up to 2-3 % are reported following CO2 dissolution into pure water at 10 MPa and up to 110 °C [2], while more moderate density increases up to ~1 % occur in brine at 12-22 MPa and 45-75 °C [3]. These density changes may produce unstable density gradients within aquifers and trigger convective flow, which substantially enhances the CO2 dissolution rate beyond that of diffusion [4,5], thus reducing the ultimate solubility trapping time scale [6].
CR IP T
The connection of density-driven convection to CO2 storage was identified early in numerical [7] and theoretical [8] studies. Enhanced CO2 dissolution rates from density-driven
convection were later verified both experimentally [5,9,10] and numerically [3,4,9]. However, changes in the aqueous phase’s density depend not only on CO2 dissolution but also on the
AN US
geothermal gradient across the formation layer, i.e., density changes are both solutally and thermally induced. Preconditions for density-driven instability can be identified by comparing the Rayleigh-Darcy number, 𝑅𝑎, (Eq. 1, section 2.2), against a critical Rayleigh-Darcy
M
number, 𝑅𝑎𝐶 .
ED
Furthermore, the 𝑅𝑎 is also related to the onset time of instability, 𝑡𝑐 , (Eq. 2, section 2.3), which determines the CO2 solubility trapping time scale. The nonlinear onset of convection
PT
occurs much later than the onset of instability [11] and depends on the initial perturbation [12,13]. The 𝑡𝑐 was studied by using linear, global or non-modal stability analysis for
CE
homogenous [14,15,16] and anisotropic [4,17-21] systems, while the onset of convection was
AC
studied numerically for homogeneous [3,11,22,23] and heterogeneous [24-26] systems. Furthermore, geochemical reactions [27], aquifer inclination [28], and exchanges with the overlying capillary transition zone [29-31] were identified to impact onset time. Several dynamic regimes of the CO2 dissolution process were identified [12,14,22,32-38]. For high 𝑅𝑎, a steady convective flux regime exists [3,25,32,33,38-40] that occurs after the onset of instability and before layer saturation. The steady convective flux, 𝐹, which governs the enhanced CO2 dissolution rate, was quantified both experimentally [32,33,37] and
ACCEPTED MANUSCRIPT numerically [25,30,38,39]. However, discrepancies exist between the obtained scaling laws of 𝐹, particularly with regards to the dependence on 𝑅𝑎. While some studies indicate the nonlinear scaling of 𝐹 with 𝑅𝑎 [32,33], other studies indicate 𝐹 to be 𝑅𝑎-independent [25,3840] but followed by a 𝑅𝑎-dependent shut-down regime [38], thus making the duration of the convective regime 𝑅𝑎-dependent. In addition, 𝐹 was found to increase when exchanges with
CR IP T
the overlying capillary transition zone are considered [30,31]. Other factors that were found to influence density-driven convection are geochemical reactions [34,41,42] and geological heterogeneity [26,43-46]. Heterogeneity affects the flow pattern, changing the flow regime and transport mechanism from convective fingering to dispersive flux with increasing
AN US
heterogeneity [26,43,44]. Thus, the onset of instability is a necessary but insufficient
condition for density-driven convection to occur. However, identifying the field conditions that are the prerequisites for density-driven instability is still paramount.
M
The thermodynamic and salinity conditions control CO2 solubility [47-49] and the
ED
subsequently induced density change, which together with the formation properties, e.g., layer thickness and permeability, determine if the prerequisites for the onset of instability exist. In
PT
addition, increases in salinity were found to delay the onset of instability [34,50], while increases in pressure and decreases in temperature cause destabilization [50]. In earlier
CE
studies, prerequisites for density-driven instability were identified for several sites
AC
[6,8,15,22,51], but the prerequisites for the onset of density-driven instability during CO2 storage have not yet been fully and systematically mapped over the broad range of conditions that are encountered in the field. This study focuses on such systematic trend analyses regarding how 𝑅𝑎, 𝑡𝑐 and 𝐹 evolve over the temperature, pressure, salinity and permeability ranges that are encountered in deep saline aquifers. Additionally, the effects of thermodynamic and salinity conditions on the convective regime duration are investigated.
ACCEPTED MANUSCRIPT Furthermore, the thermal contribution to convection, which is induced by the geothermal gradient, was previously found to be negligible compared to the solutal contribution [8,46,51]. However, we extend the analysis in this study to include both direct and indirect geothermal effects while analysing the influences of the geothermal gradient (i.e., basin type) on the trends in the onset of instability, 𝑡𝑐 and 𝐹.
CR IP T
In addition, the relative influence of different field conditions and properties - such as storage depth, basin type, porosity, permeability, layer thickness and salinity - on density-driven
instability and convective mixing is identified through local and global sensitivity analyses.
AN US
The theoretical findings are discussed in relation to site-specific cases in the form of densitydriven instability prerequisites for CO2 pilot test and storage sites.
For this study, a new analysis tool was developed with the aim to avoid using simplifying
M
assumptions that are commonly made during the estimation of 𝑅𝑎. The tool was used to rigorously calculate 𝑅𝑎, 𝑡𝑐 and 𝐹 by fully accounting for the temperature, pressure and
ED
salinity dependency of CO2 solubility, induced density changes, diffusivities and viscosity.
PT
This approach was possible through the implementation of state-of-the-art solubility, density and viscosity models and algorithms for calculating diffusivities.
CE
These findings give new insight into the field conditions that favour or hinder enhanced
AC
solubility trapping and identify the key field parameters that affect the process; thus, they are expected to be useful during the selection of appropriate CO2 storage sites and for capacity estimation and risk assessment. In addition, a simple and efficient procedure to obtain an accurate estimate of site-specific solutal 𝑅𝑎 is presented.
2. Methods 2.1 Description of the system
ACCEPTED MANUSCRIPT The conceptual model consists of an isotropic, homogenous porous medium of thickness 𝐻 (Fig. 1). Injected CO2 has migrated to the upper part of the system because of buoyancy. Two regions are distinguished: an upper multi-phase region that contains CO2-rich phase and CO2saturated aqueous phase, and an underlying single-phase region that contains aqueous phase with dissolved CO2. The CO2 that is dissolved in the single-phase region induces density
CR IP T
changes in the aqueous phase and, potentially, instability. Several upper boundary conditions are considered: impermeable, partially or fully permeable to flow. These correspond to
different amounts of exchange with the capillary transition zone. In the single-phase region, the CO2 concentration is initially zero (although in some studies assumed to vary linearly with
AN US
depth, see Table 1), while the temperature profile is linear and set by the geothermal gradient.
2.2 Onset of Instability
M
Prerequisites for density-driven instability exist when 𝑅𝑎 > 𝑅𝑎𝐶 . 𝑅𝑎, which is given by Eq. 1,
number, 𝑅𝑎𝑆 [51,52]:
𝑤ℎ𝑒𝑟𝑒
PT
𝑅𝑎 = 𝑅𝑎 𝑇 + 𝑅𝑎𝑆
ED
is the sum of the thermal Rayleigh-Darcy number, 𝑅𝑎 𝑇 , and the solutal Rayleigh-Darcy
𝑅𝑎 𝑇 =
𝑔𝑘∆𝜌𝑇 𝐻 , 𝜇𝜅
𝑅𝑎𝑆 =
𝑔𝑘∆𝜌𝑆 𝐻 𝜇𝜑𝐷
(1)
CE
Here, 𝑔 is the magnitude of gravitational acceleration, 𝑘 is the permeability, 𝐻 is the layer thickness, 𝜇 is the formation fluid viscosity, 𝜑 is the porosity, and ∆𝜌𝑆 and ∆𝜌𝑇 are,
AC
respectively, the solutally and thermally induced density changes. 𝐷 and 𝜅 represent the molecular and thermal diffusivity, respectively. Note that 𝜇, ∆𝜌𝑆 , ∆𝜌𝑇 , 𝐷 and 𝜅 are all temperature, pressure and salinity dependent. Values of 𝑅𝑎𝐶 for different boundary conditions of the system in Fig. 1 are given in Table 1. 𝑅𝑎𝐶 decreases with increasing permeability of the upper boundary. In this study, the 𝑅𝑎𝐶 range 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 was adopted. This 𝑅𝑎𝐶
ACCEPTED MANUSCRIPT range, although 𝑅𝑎𝐶 nearly triples, can be considered narrow when compared to the 𝑅𝑎 range that was studied.
The linear onset time of instability, 𝑡𝑐 , is given by Eq. 2 [4]:
𝑡𝑐 = 𝑐0
1 𝐻2 (𝜇𝜑)2 𝐷 ∙ = 𝑐 0 (𝑅𝑎𝑆 )2 𝐷 (𝑔𝑘∆𝜌𝑆 )2
CR IP T
2.3 Onset time of Instability
(2)
AN US
The value of the numerically derived constant 𝑐0 was reported [3,12,14,15,17-
19,22,23,25,29,30] to range from 𝑐0 =31 [30] to 𝑐0 =5000 [25], depending on the method, initial conditions and stability criteria that were used [22]. In this study, the lower and upper limits of 𝑐0 =31 and 𝑐0 =146, which correspond to full exchange or no exchange with the
ED
2.4 Steady convective flux
M
capillary transition zone [30], were used.
PT
Scaling laws for the non-dimensional steady convective flux, 𝐹 ∗ , are given in Table 2. Both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws are included. Values of 𝐹 ∗ for several
CE
upper boundary conditions with different permeabilities are given. The dimensional steady
AC
convective flux, 𝐹, is given by Eq. 3 [37]: 𝐹 = 𝐹∗ ∙
𝑔𝑘∆𝜌𝑆 ∙ 𝐶𝑠𝑎𝑡 𝜇
(3)
Here, 𝐶𝑠𝑎𝑡 is the mass of solute per unit volume of saturated solution. The onset time of convection depends on the initial perturbation size [12]. A numerically based parameterization that approximates the steady convective regime and its duration was recently presented by Slim [38] and is included in Table 2.
ACCEPTED MANUSCRIPT
2.5 Developed analysis tool A comprehensive 𝑅𝑎-𝑡𝐶 -𝐹-calculation model was developed and coded in MATLAB. See Fig. 2 for the model structure with submodels and work flow. To ensure the accurate and continuous calculation of 𝑅𝑎, 𝑡𝐶 , and 𝐹, the analysis tool incorporates submodels that account
CR IP T
for the temperature, pressure and salinity dependence of CO2 solubility into brine, thermally and solutally induced density changes, the viscosity of the brine and the thermal and
molecular diffusivity. The analysis tool is valid for a broad range of temperatures, 30-130 °C; pressures, 7.85-60 MPa; and salinities, 0-6 m NaCl (mol NaCl/kg H2O). This salinity range,
AN US
which is equivalent to 0-350,000 ppm, was previously applied by Kumar et al. [53], who studied geological CO2 storage trapping mechanisms. 2.5.1 Solubility submodels
M
Two solubility models, Duan and Sun [47] and Spycher and Pruess [48,49], were incorporated
ED
into the analysis tool. Both models are widely used to calculate CO2 solubility under CO2 storage conditions. One of the solubility models, [49], is applicable for temperatures of 12-
PT
100 °C, pressures of 0-60 MPa and salinities of 0-6 m NaCl, while the other model, [47], is valid for temperatures of 0-260 °C, pressures of 0-200 MPa and salinities of 0-4.3 m NaCl.
CE
The use of both solubility models provides a comparison opportunity. The CO2 solubility, 𝑥𝑆 ,
AC
increases with pressure and decreases with ionic strength, while the temperature dependence is non-monotonic. At low pressures, such as <10 MPa, solubility decreases with temperature. At higher pressures, solubility decreases with temperature until a pressure dependent minimum and then increases with temperature [47]. 2.5.2 Density submodel
ACCEPTED MANUSCRIPT A slightly modified version of the aqueous phase density model by Li et al. [54], which accounts for dissolved CO2, was incorporated into the analysis tool. The model by Li et al. [54] combines the density model by Mao and Duan [55] for aqueous NaCl solutions with a similar approach to that of Duan et al. [56] for the density of CO2-H2O-NaCl solutions. Our modifications consisted of adopting the thermodynamic formulation IAPWS97 [57] for pure
CR IP T
water density and calculating the volumetric Debye-Hückel limiting law slope by the linearization of an extensive dataset [58] of tabulated values. Furthermore, the analysis tool calculates ∆𝜌𝑆 as the density difference between the CO2-H2O-NaCl solution and the resident brine. The ∆𝜌𝑇 is calculated as the density difference between resident brine at the storage
AN US
layer top and storage layer bottom. The latter is determined by using, e.g., the formation thickness, 𝐻, and the geothermal gradient, 𝛻𝑇𝑔𝑒𝑜 , to derive the bottom temperature and pressure.
M
2.5.3 Viscosity submodel
ED
The dynamic viscosity model by Mao and Duan [59] was incorporated into the analysis tool to calculate the viscosity of NaCl solutions as a function of temperature, pressure and salinity.
PT
The viscosity model was chosen because of its applicability to the pressure and temperature
CE
conditions that are found under geological CO2 storage and a wide salinity range of 0-6 m NaCl. Viscosity is highly temperature dependent, decreasing with increasing temperature.
AC
Salinity also affects brine viscosity but to a lesser extent than temperature. The salinity effect on viscosity depends on the electrolyte; for NaCl solutions, the viscosity increases with ionic strength [59]. The influence of pressure is comparably small, as higher pressure raises the brine viscosity only slightly [60]. This study only addresses density-driven instability that is induced by ∆ρS and ∆ρT , and potential influences on instability that arise from viscosity changes that are caused by dissolved CO2 are not considered. Viscosity changes that follow CO2 dissolution are assumed to be negligible.
ACCEPTED MANUSCRIPT 2.5.4 Molecular diffusivity submodel The CO2 diffusivity, 𝐷, flattens concentration gradients within the porous medium. In this study, a dilute solution was assumed and the Wilke-Chang correlation [61], where the molecular diffusivity changes with temperature and solvent viscosity, 𝜇, was applied, as in
CR IP T
Eq. 4: ̃𝑤 )1/2 𝑇 7.4 ∙ 10−8 (𝜃𝑀 𝐷= 𝜇𝑉̅𝑆0.6
(4)
̃𝑤 is Here, 𝜃 is an empirical constant of the solvent, in this case water, with a value of 2.6. 𝑀
AN US
the molecular weight of water, and T is the temperature in Kelvin. The molar volume of the solute, i.e., CO2, 𝑉̅𝑆 , is 34.0 cm3/mol [61,62]. 2.5.5 Thermal diffusivity submodel
M
The thermal diffusivity, 𝜅, smooths temperature differences. 𝜅 is primarily temperature
ED
dependent, decreasing with temperature. The 𝜅 that is used to calculate 𝑅𝑎 is defined as the ratio between the thermal conductivity of the medium, 𝜆, and the fluid heat capacity, (𝑐𝜌)𝑓 ,
𝜅 = 𝜆 / (𝑐𝜌)𝑓
(5)
CE
PT
as in Eq. 5 [52]:
AC
Here, 𝑐𝑓 and 𝜌𝑓 are the heat capacity and density of the fluid, respectively. In this study, the heat capacity and density of water were calculated by employing the thermodynamic formulation of IAPWS97 [57], taking into account both pressure and temperature. The thermal conductivity, 𝜆, was estimated with the correction [63,64] in Eq. 6: 𝜆(𝑇) = 𝜆20°𝐶 − 10−3 ∙ (𝑇 − 293) ∙ (𝜆20°𝐶 − 1.38) ∙ [𝜆20°𝐶 ∙ (1.8 ∙ 10−3 𝑇)−0.25∙𝜆20°𝐶 + 1.28] ∙ 𝜆20°𝐶 −0.64
(6)
ACCEPTED MANUSCRIPT Here, the temperature, T, is given in Kelvin and the thermal conductivity at 20 °C, 𝜆20°𝐶 , is given in W/(m K). A sedimentary basin of sandstone with 𝜆20°𝐶 = 1.5 W/(m K) was assumed, a value in the range that was reported for sedimentary rocks [65]. Eq. 6 was found to be valid for both dry and liquid-saturated sandstones with porosities that are not too high [63]. In this study, liquid-saturated sandstone was considered.
CR IP T
2.5.6 Calculation of 𝑅𝑎, 𝑡𝐶 and 𝐹
The values of 𝜇, ∆𝜌𝑆 , ∆𝜌𝑇 , 𝐷 and 𝜅, which were calculated by the above submodels, were transferred (Fig. 2) to submodels to compute 𝑅𝑎, 𝑅𝑎𝑆 and 𝑅𝑎 𝑇 (Eq. 1), 𝑡𝐶 (Eq. 2) and 𝐹
AN US
(Table 2 and Eq. 3). All the scaling laws that are given in Table 2 are incorporated into the analysis tool.
2.6 Comparison study
M
The performance of each submodel that was described in sections 2.5.1 through 2.5.5 was
ED
verified against its original source study with excellent agreement. Furthermore, the analysis tool’s performance was compared to the study by Lindeberg and Wessel-Berg [8] as a
PT
benchmark by comparing calculated S-numbers for four different North Sea reservoir
CE
conditions. The S-number definition is given in the figure caption of Fig. 3. Because Lindeberg and Wessel-Berg [8] assumed a constant ∆𝜌𝑆 , a first comparison was made by
AC
replacing the temperature-, pressure- and salinity-dependent ∆𝜌𝑆 that was calculated by the analysis tool with their constant value. Fig. 3 shows that the calculated S-numbers concurred with the S-numbers that were estimated by Lindeberg and Wessel-Berg [8]. Only slight differences originated from the marginal differences in the thermal and molecular diffusivities and viscosities and the differences in the estimated ∆ρT . In the present study, lower values of ∆ρT were obtained as the analysis tool also accounts for the effect of the pressure gradient, which counteracts thermal expansion.
ACCEPTED MANUSCRIPT
To evaluate the effect of a non-constant ∆𝜌𝑆 -assumption, the pressure-, temperature- and salinity-dependent ∆𝜌𝑆 was calculated by the analysis tool and compared to the constant ∆𝜌𝑆 . The S-numbers that were obtained by the analysis tool showed similar trends as those that were obtained by Lindeberg and Wessel-Berg [8] but were systematically lower. This
CR IP T
indicates that the constant ∆𝜌𝑆 -assumption by Lindeberg and Wessel-Berg [8] overestimates S-numbers (case 1 by 21.5 %, case 2 by 40.1 %, case 3 by 52.2 % and case 4 by 57.4 %) for all conditions. The largest overestimations coincide with the highest temperature and pressure
AN US
conditions.
2.7 Sensitivity analysis
For this study, local and global sensitivity analyses were applied to identify influential
M
parameters on system variables. The analyses were conducted by connecting the developed
ED
analysis tool to iTOUGH2 [66] through the PEST protocol [68,69]. In the local sensitivity analysis, partial derivatives, i.e., the ratios between variable changes
PT
and the parameter perturbations that give rise to them, were scaled to obtain dimensionless
CE
sensitivity coefficients. The scaling was conducted by multiplying the partial derivatives with the ratio of the expected parameter variation and variable standard deviation. The
AC
dimensionless sensitivity coefficients permitted a comparison of the relative influence of parameters of different types, units or scales [66]. In addition, a global sensitivity analysis was conducted to investigate the parameter space, calculate sensitivity measures and identify parameter-interaction effects. The Morris one-at-a-time (OAT) method [67] was used, which partitions the parameter ranges into equally sized intervals and randomly generates sets of parameters. The generated parameter values are then perturbed and the changes in the variables are calculated. The analysis is performed over multiple “paths” (sets of parameter
ACCEPTED MANUSCRIPT values and orders of perturbing parameters) to obtain a global sensitivity index, the mean elementary effect (EE), for each parameter over the parameter space. A high absolute EE indicates a highly influential parameter, while a high standard deviation (SD) of the EE indicates the presence of non-linearity or interaction effects [66]. Furthermore, the uncertainty
CR IP T
in an EE estimate is given by 𝑆𝐷/𝑟 0.5, where 𝑟 is the number of paths [67].
2.8 Investigated cases
Analyses of the prerequisites for density-driven instability were conducted over broad ranges of conditions that are likely found under CO2 storage in deep saline aquifers, such as
AN US
temperatures of 30-130 °C, pressures of 7.85-60 MPa and aqueous NaCl solutions of 0-6 m (mol NaCl/kg H2O).
Trends in the CO2 solubility and subsequent ∆𝜌𝑆 over these broad ranges of conditions were
M
studied by using the analysis tool. Then, two aspects were studied to obtain an overview of
ED
the relative solutal and thermal contributions to density-driven instability over the ranges of these conditions. First, the trends of the relative solutal and thermal driving forces of
PT
convection, i.e., the buoyancy ratio, 𝑁 = ∆𝜌𝑆 /∆𝜌𝑇 , were studied within aquifers with different geothermal gradients, salinities and aquifer properties (Table 3). Second, the relative
CE
solutal and thermal contributions to the onset of instability, i.e., 𝑅𝑎, were studied through trends in 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 for the ranges of conditions that are given in Table 3. The ratio 𝜑𝑅𝑎𝑆 /
AC
𝑅𝑎 𝑇 is influenced by both 𝑁 and the diffusivity ratio. Thereafter, systematic calculations of 𝑅𝑎 over the broad thermodynamic and salinity ranges and permeabilities, which are given in Table 3, were conducted to obtain an overview of the 𝑅𝑎 range that is expected during CO2 storage in deep saline aquifers and to investigate the 𝑅𝑎 trends and limiting conditions for density-driven instability.
ACCEPTED MANUSCRIPT Furthermore, the impact of basin type on the prerequisites for density-driven instability was investigated; the 𝑅𝑎 and corresponding 𝑡𝐶 were calculated and analysed for the thermodynamic and salinity conditions that are likely found in a cold basin and a warm basin, respectively. Cold basins have low surface temperatures and geothermal gradients compared to warm basins, thus generally exhibiting lower temperatures at the same depths [70]. In
CR IP T
earlier studies [70,71], the basin type was shown to affect the CO2 storage capacity of free phase CO2; however, its impact on density-driven convective mixing, the focus of this study, has not been explored. The values of the surface temperature and geothermal gradients that are used in the present study (Table 3) were obtained from Bachu [71] and the 𝑡𝐶 by using Eq.
AN US
2.
To relate this study's theoretical findings to real cases, the 𝑅𝑎𝑆 and its corresponding 𝑡𝐶 for five CO2 pilot test and storage sites were calculated by using the parameters that are shown in
M
Table 4 and Eq. 2.
ED
Additionally, a procedure that permits the application of this study's results in a straightforward manner was developed to estimate 𝑅𝑎𝑆 for sites of one's own choosing while
PT
accounting for salinity and thermodynamic conditions. To create a basis for accurate and
CE
efficient site-specific estimations of 𝑅𝑎𝑆 , the thermodynamic- and salinity-dependent part of Eq. 1, i.e., ∆𝜌𝑆 /𝜇𝐷, hereafter referred to as the dynamic factor, was generated for a wide
AC
range of temperature (30-130 ˚C), pressure (10-60 MPa) and salinity (0-4 m NaCl) conditions. The impacts of the thermodynamic and salinity conditions and the basin type on the trends in 𝐹 and the duration of the convective regime were also studied. A comparison of 𝐹 that used both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws under cold basin conditions, which are given in Table 3, was conducted. The value of k that was used made the comparison of the scaling laws possible within their valid 𝑅𝑎𝑆 ranges, shown in Table 2. For the comparison, the
ACCEPTED MANUSCRIPT salinity was set to 0 m NaCl, which led to upper limit estimates of 𝐹. The impacts of the basin type and storage depth on the trends in 𝐹 were studied for three different salinities by using the scaling law by Slim [38], Eq. 3 and the properties in Table 3. The impacts of the basin type, salinity and storage depth on the duration of the convective regime were also
CR IP T
investigated by using the parameterization by Slim [38] in Table 2. Additionally, local and global sensitivity analyses were conducted to determine the relative influence of different field conditions and properties on density-driven instability and
convective mixing. In the local sensitivity analysis, the influences of the field parameters
AN US
(storage depth, basin type, i.e., ∇𝑇𝑔𝑒𝑜 , 𝜑, 𝑘, 𝐻 and salinity) on the variables 𝑥𝑆 , ∆𝜌𝑆 , 𝐷, 𝜇, 𝑅𝑎𝑆 , 𝐹, and 𝑡𝐶 and the duration of the convective regime were investigated. The analysis included 𝐹 that was based on both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws. However, only the influences of the field parameters on 𝑅𝑎𝑆 were analysed in the global
ED
analyses are given in Table 5.
M
sensitivity analysis (Morris OAT method). The parameters that were used in the sensitivity
PT
3. Results
CE
3.1 Trends in CO2 solubility and ∆𝜌𝑆 under broad ranges of conditions The CO2 solubility into brine and the subsequent ∆𝜌𝑆 % at different temperatures, pressures
AC
and salinities are shown in Fig. 4 and Fig. 5, respectively. Similar trends in CO2 solubility and ∆𝜌𝑆 % are produced (although the results deviate slightly at high temperatures and pressures at low salinities) by both solubility models [47,49]: decreasing trends with decreasing pressure and increasing salinity. Furthermore, the CO2 solubility first decreases with temperature until a pressure-dependent minimum is reached and then increases with temperature, while ∆𝜌𝑆 % decreases with increasing temperature. However, the impact of the temperature and pressure
ACCEPTED MANUSCRIPT conditions on both CO2 solubility and ∆𝜌𝑆 % diminishes with increasing salinity. Because of the salinity and temperature validity limits of the solubility models, only results from the Spycher and Pruess [49] or Duan and Sun [47] solubility model are shown for salinities over 4 m NaCl or temperatures over 100 °C, respectively. Fig. 5 shows that ∆𝜌𝑆 % is generally around 1 % to a few tenths of a percentage for the wide range of conditions that were studied.
CR IP T
The maximum ∆𝜌𝑆 %, ~2 % (not shown), corresponds to thermodynamic conditions of 30 °C, 60 MPa and 0 m NaCl. However, at high temperatures and salinities, ∆𝜌𝑆 % is negative, i.e., a density decrease follows CO2 dissolution into brine. This result is consistent with the result by Lu et al. [74] on CO2 storage efficiency, who showed the temperature, or “equal density
AN US
temperature”, beyond which CO2-saturated brine is lighter than resident brine (∆𝜌𝑆 < 0) to be 118 °C at 20 MPa and 4 m salinity. However, for the wider salinity range up to 6 m NaCl that was investigated in the present study, Fig. 5 shows that the “equal density temperature”
M
decreases with increasing salinity for a certain pressure; at conditions of 20 MPa and 6 m
ED
NaCl, the “equal-density temperature” is 80 °C and is even lower at lower pressures.
instability
PT
3.2 The relative solutal and direct thermal contributions to the onset of
CE
Previous studies [8,46,51] have found that the thermal contribution to convection is negligible compared to the solutal contribution; ∆𝜌𝑇 is generally one order of magnitude less than ∆𝜌𝑆
AC
[51] and 𝑅𝑎 𝑇 is orders of magnitude less than 𝑅𝑎𝑆 [8]. In the present study, a relative comparison of ∆𝜌𝑆 and ∆𝜌𝑇 , i.e., 𝑁, was made within an aquifer for different ∇𝑇𝑔𝑒𝑜 and salinities (Fig. 6). The assumed 𝐻 in this study (100 m) is within the range of 𝐻 that is considered in previous studies [8,51]. The 𝑁 varies widely, depending on the thermodynamic and salinity conditions and increases as the ∇𝑇𝑔𝑒𝑜 and salinity decrease. The highest ratios are obtained for thermodynamic conditions that correspond to shallow depths. Over the wide
ACCEPTED MANUSCRIPT ranges of conditions that were studied, ∆𝜌𝑇 varies from one order of magnitude smaller than ∆𝜌𝑆 to reaching equal magnitude or even larger. At larger depths, the ratios, which vary with ∇𝑇𝑔𝑒𝑜 , turn negative as a result of the above-mentioned induced density decrease. Fig. 7 shows 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 for the same ∇𝑇𝑔𝑒𝑜 and salinities that were used to study 𝑁. The trend in φRaS /RaT with depth is a combined effect of the trends in both N and the diffusivity ratio.
CR IP T
These trends are similar, but the influence of the latter is greater. 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 increases with decreasing ∇𝑇𝑔𝑒𝑜 , depth, 𝜑 and salinity. The ratio is generally large, with 𝑅𝑎𝑆 being several hundreds to tens of thousand times larger than 𝑅𝑎 𝑇 . However, the ratio turns negative at a
AN US
certain depth, which is determined by ∇𝑇𝑔𝑒𝑜 , and high salinities.
3.3 Trends in 𝑅𝑎 under broad ranges of conditions
In Fig. 8, 𝑅𝑎 values for different temperature, pressure, salinity and permeability conditions
M
are shown. The Y-axis scale of each graph is different. The displayed trends in 𝑅𝑎 are of greater importance than the absolute values because the absolute values are based on likely,
ED
but assumed, formation properties such as 𝐻 and 𝜑. The trends in 𝑅𝑎 are similar to those that
PT
are displayed by ∆𝜌𝑆 % in Fig. 5: a decreasing trend with increasing temperature and salinity and an increasing trend with increasing pressure. The pressure influence on 𝑅𝑎 diminishes
CE
with temperature at high salinities (≥4 m NaCl). Negative values of 𝑅𝑎 correspond to negative
AC
∆𝜌𝑆 %, which is mentioned in connection to Fig. 5. The solid red lines in Fig. 8 show the 𝑅𝑎𝐶 range 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 . Values of 𝑅𝑎 that fall above the upper limit of the 𝑅𝑎𝐶 range indicate prerequisites for the onset of instability under the specific thermodynamic and salinity conditions. The 𝑘 has a large influence on the magnitude of 𝑅𝑎. In relatively high 𝑘 (10-12 m2) and low-moderate salinity aquifers, 𝑅𝑎 is generally large, on the order of a thousand to tens of thousands, i.e., 𝑅𝑎 ≫ 𝑅𝑎𝐶 , and the prerequisites for the onset of densitydriven instability are good. For even higher 𝑘 or thicker layers, 𝑅𝑎 values of a few hundreds
ACCEPTED MANUSCRIPT of thousands can be expected. However, even in aquifers with high 𝑘, 𝑅𝑎 values that are close to zero or negative can be found at high salinity and certain thermodynamic conditions. For aquifers with lower k (10-14 m2) and low-moderate salinity, 𝑅𝑎 is one order higher or of the same order of magnitude as the RaC range. Under these conditions, instability is highly dependent on the prevailing thermodynamic and salinity conditions. In low 𝑘 (10-15 m2)
CR IP T
aquifers, the 𝑅𝑎 values are generally small, falling below or within the 𝑅𝑎𝐶 range, and the prerequisites for the onset of density-driven instability are limited. Furthermore, Fig. 8 shows that 𝑅𝑎 values in high salinity aquifers are generally low, except at high 𝑘 and low-moderate temperatures. In addition, large negative values of 𝑅𝑎 are obtained for high salinity and
AN US
temperature conditions at high 𝑘.
3.5 Trends in 𝑅𝑎 and 𝑡𝐶 under field conditions - impact of the basin type (i.e.,
M
∇𝑇𝑔𝑒𝑜 )
ED
In Fig. 9, 𝑅𝑎 values for different 𝑘 and salinities for a cold basin at 800-3500 m depth and a warm basin at 800-2400 m depth are shown. For the warm basin case, results were obtained
PT
down to a depth of 1800 m by using the Spycher solubility model [49] and at lower depths down to 2400 m by using the Duan and Sun solubility model [47] because of the temperature
CE
validity limits of the solubility models that were used. Furthermore, negative 𝑅𝑎 values are
AC
not shown because of the logarithmic Y-axis scales. A general trend is that the highest 𝑅𝑎 is obtained for temperature and pressure conditions that correspond to shallow depths of around 800 m, and then 𝑅𝑎 decreases with depth. Furthermore, the 𝑅𝑎 values are generally higher in the cold basin than in the warm basin because of the smaller ∇𝑇𝑔𝑒𝑜 of the former, which induces lower temperatures and thus larger ∆𝜌𝑆 (Fig. 5) and leads to larger 𝑅𝑎 in turn. In the cold basin, 𝑅𝑎 generally exceeds the
ACCEPTED MANUSCRIPT assumed 𝑅𝑎𝐶 range at high 𝑘 (10-13 m2 to 10-12 m2), indicating prerequisites for the onset of instability, except at high salinities. At high salinities, the depth determines the prerequisites for density-driven instability. Similar trends are seen in the warm basin, but the depths at which the prerequisites for instability exist are shallower (at even lower salinities) than in the cold basin. For aquifers in the intermediate 𝑘 range (10-14 m2 to 10-13 m2) with low salinities
CR IP T
(≤1 m NaCl) and for both basin types, the 𝑅𝑎 values generally exceed the 𝑅𝑎𝐶 range for the studied depths. However, as the salinity increases, the influence of depth on the onset of
instability increases. In lower 𝑘 (10-15 m2 to 10-14 m2) aquifers, the prerequisites for densitydriven instability are more limited and are highly affected by the salinity, depth and basin type
AN US
(i.e., ∇𝑇𝑔𝑒𝑜 ). For formations with low 𝑘 (≤ 10-15 m2), the 𝑅𝑎 values are small. Only for low salinity and shallow depths in the cold basin do the 𝑅𝑎 values fall within the considered 𝑅𝑎𝐶 range but do not exceed it, i.e., the prerequisites for instability are very limited for low 𝑘.
M
In Fig. 10, the impact of salinity on 𝑅𝑎𝑆 is shown under cold and warm basin conditions. The
ED
𝑅𝑎𝑆 values that are obtained under cold basin conditions generally exceed the 𝑅𝑎𝐶 range for low-moderate salinities. However, for high salinities, the depth, i.e., thermodynamic
PT
conditions, highly influences whether prerequisites for instability exist. For a cold basin with a salinity of ~5 m NaCl, the temperature and pressure conditions for which 𝑅𝑎𝑆 first falls
CE
within the 𝑅𝑎𝐶 range correspond to a depth of ~3300 m but decreases to 2600 m for a salinity
AC
of ~6 m NaCl. Under warm basin conditions, the depths at which 𝑅𝑎𝑆 first falls within the 𝑅𝑎𝐶 range are even shallower than those under cold basin conditions and occur at even lower salinities. In the warm basin, only salinities <3 m NaCl give rise to values of 𝑅𝑎𝑆 that exceed the 𝑅𝑎𝐶 range within the whole depth range. As in the cold basin, the depth has a major influence on the prerequisites for the onset of instability at higher salinities. For the warm basin and salinities of ~4 m NaCl, ~5 m NaCl and ~6 m NaCl, 𝑅𝑎𝑆 first falls within the 𝑅𝑎𝐶 range at a depth of approximately 1800 m, 1400 m and 1100 m, respectively.
ACCEPTED MANUSCRIPT In Fig. 11, the 𝑡𝐶 values for different salinities under cold and warm basin conditions are shown. These values of 𝑡𝐶 are based on the values of 𝑅𝑎𝑆 in Fig. 10. Only 𝑡𝐶 for 𝑅𝑎𝑆 ≥13.8 are shown. Furthermore, the values of 𝑡𝐶 were calculated by using Eq. 2 and the numerical constant 𝑐0 =31, which represents impermeable upper boundary conditions [30], thus giving the minimum estimates of 𝑡𝐶 compared to assuming a permeable upper boundary condition.
CR IP T
One trend (Fig. 11) is that 𝑡𝐶 increases with salinity and depth and the influence of the depth increases with salinity. However, 𝑡𝐶 is generally smaller for cold basin conditions than for warm basin conditions at the same depth.
AN US
3.6 Field sites and a new procedure for the accurate site-specific estimation of 𝑅𝑎𝑆
The calculated values of 𝑅𝑎𝑆 and 𝑡𝐶 for five CO2 pilot test and storage sites, both off- and
M
onshore, are given in Table 6. These pilot test and storage sites have different temperature, pressure and salinity settings and formation properties, as given in Table 4. The largest 𝑅𝑎𝑆 is
ED
obtained for the Sleipner site; the ∆𝜌𝑆 % is almost 1 %, considerably higher than for the other
PT
sites. The 𝑅𝑎𝑆 for Frio, Ketzin and Snøhvit are of the same order of magnitude. The lowest 𝑅𝑎𝑆 is obtained for the In Salah site in the Krechba field. This 𝑅𝑎𝑆 , which falls below the 𝑅𝑎𝐶
CE
range, occurred because of the combined effect of small values of 𝑘, 𝐻, and ∆𝜌𝑆 under the thermodynamic and salinity conditions of the site. The 𝑡𝐶 values for the sites are given in
AC
Table 6. The difference in the 𝑡𝐶 between the Ketzin and Snøhvit sites, despite the sites’ similar 𝑅𝑎𝑆 and 𝐻 values, arises from the difference in 𝐷. A two-step procedure was developed to apply this study's results in a straight-forward manner to estimate 𝑅𝑎𝑆 for a site of one's own choosing: 1) Given the temperature, pressure and salinity of a specific site, Fig. 12 can be used to obtain the dynamic factor, i.e., ∆𝜌𝑆 /𝜇𝐷 in Eq. 1, for a wide range of conditions (proxy-equations for the dynamic factor can be found in
ACCEPTED MANUSCRIPT Appendix A), and approximations of the dynamic factor for salinities that are not shown can be made by interpolation between two dynamic factors that are obtained for a lower and higher salinity. 2) Once the dynamic factor is obtained, the site-specific 𝑅𝑎𝑆 (Eq. 1) can be estimated by dividing by 𝜑 and multiplying with 𝐻 (in unit m), 𝑘 (in unit m2) and 𝑔 (in unit m/s2), which varies with layer inclination. The estimated 𝑅𝑎𝑆 can then serve as an indicator
CR IP T
for the prerequisite for the onset of density-driven instability at a site and used to evaluate the suitability of storage sites.
3.7 Trends in 𝐹 under field conditions and impacts on the duration of the
AN US
convective regime
In Fig. 13, the 𝐹 values that are predicted by different scaling law types (𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent) and assumed upper boundary conditions (Table 2) are shown under cold
M
basin conditions. Both the non-dimensional mass flux, 𝐹 ∗ , and corresponding dimensional mass flux, 𝐹, are shown. The mass fluxes are obtained for zero salinity conditions and
ED
therefore constitute an upper bound estimate of the flux for the given thermodynamic
PT
conditions. Fig. 13 gives an overview of the implications of using a specific type of scaling law or upper boundary condition on the predicted 𝐹. The comparison shows consistency in
CE
the predicted trend in 𝐹 whether 𝑅𝑎𝑆 -dependent or 𝑅𝑎𝑆 -independent scaling laws were used; 𝐹 increased with increasing depth. In Fig. 13, the experimentally derived 𝑅𝑎𝑆 -dependent
AC
scaling laws are shown with their uncertainty intervals (dotted lines). For the range of 𝑅𝑎𝑆 values that were studied, the steady convective flux that was predicted by the 𝑅𝑎𝑆 independent scaling laws for impermeable (by, e.g., Slim [38] and Elenius et al. [30]) and partially (0.2) permeable upper boundary fall within the uncertainty interval of the 𝑅𝑎𝑆 dependent scaling law that was experimentally derived by Neufeld et al. [32]. A relatively wide range of predicted flux sizes is shown in Fig. 13 for different upper permeability
ACCEPTED MANUSCRIPT boundary conditions. The smallest steady convective fluxes result from the 𝑅𝑎𝑆 -dependent scaling law that was experimentally derived by Backhaus et al. [33]. Results are only shown for cold basin conditions because similar results were obtained for warm basin conditions. The effects of salinity, basin type and storage depth on the trends in 𝐹 (using the scaling law by Slim [38]) down to a depth of 2500 m and the corresponding ∆𝜌𝑆 /𝜇 ratio are shown in Fig.
CR IP T
14. Only values of 𝐹 at conditions for which the system’s 𝑅𝑎𝑆 is within the valid 𝑅𝑎𝑆 range (Table 2) are shown in Fig. 14, which is the reason why 𝐹 under warm basin conditions (at 3 m NaCl) is not shown for the entire depth range. A substantial decrease in 𝐹 is visible as
AN US
salinity increases in both cold and warm basins, but the basin type (i.e., ∇𝑇𝑔𝑒𝑜 ) also has a visible impact on 𝐹. Furthermore, the trends in 𝐹 closely follow the salinity-dependent trends in the ∆𝜌𝑆 /𝜇 ratio. At low salinities (0-1 m NaCl), 𝐹 and ∆𝜌𝑆 /𝜇 increase with depth. However, as salinity increases, the trend in 𝐹 with depth also changes. The trend in 𝐹 evolves,
M
first increasing with depth, then becoming more uniform, and finally decreasing with depth as
ED
salinity increases.
In addition, Fig. 15 shows the effects of salinity, basin type and storage depth on the
PT
convective flux regime with regards to 𝐹 and the 𝑅𝑎𝑆 -dependent duration (when applying the
CE
parameterization by Slim [38]). Only the convective flux regime is shown, not the preceding or successive flux regimes. Fig. 15 shows that combinations of factors such as different
AC
salinity, storage depth and basin type (i.e., ∇𝑇𝑔𝑒𝑜 ) highly affect 𝐹 and the approximate onset time and duration of the convective regime; e.g., high salinity, large storage depth and warm basin conditions (large ∇𝑇𝑔𝑒𝑜 ), which correspond to a low 𝑅𝑎𝑆 , delay the onset of convection, decrease the magnitude of the steady convective flux and shorten the duration of the convective regime.
3.8 The relative influence of different field conditions and properties
ACCEPTED MANUSCRIPT The local sensitivity analysis, where scaled sensitivity coefficients were used to compare the influence of different field conditions and properties on the 𝑥𝑆 , ∆𝜌𝑆 , 𝐷, 𝜇, 𝑅𝑎𝑆 , 𝐹, 𝑡𝐶 and duration of the convective regime, has provided several findings (Fig. 16). The analysis indicates that 1) 𝑅𝑎𝑆 is sensitive to all the studied parameters. However, 2) 𝑘, 𝜑 and 𝐻 have the highest influence on 𝑅𝑎𝑆 . 3) The 𝑥𝑆 , ∆𝜌𝑆 , 𝐷 and 𝜇 are sensitive to ∇𝑇𝑔𝑒𝑜 , salinity and
CR IP T
storage depth. 4) The most influential parameters on 𝐹 when applying both 𝑅𝑎𝑆 -independent (𝐹 a ) and 𝑅𝑎𝑆 -dependent (𝐹 b ) scaling laws are 𝑘 and salinity. However, 5) for 𝐹 a, the
sensitivity to 𝜑 and 𝐻 differs, which is expected because 𝑅𝑎𝑆 is highly sensitive to these parameters. 6) The ∇𝑇𝑔𝑒𝑜 and storage depth seem to have an influence on 𝐹. 7) The most
AN US
influential parameters on 𝑡𝐶 are 𝑘 and 𝜑, followed by salinity, ∇𝑇𝑔𝑒𝑜 and storage depth. 8) All the studied parameters influence the duration of the convective regime, but the most influential are 𝑘, 𝐻, 𝜑 and salinity.
M
The result from the global sensitivity analysis (Morris OAT method) regarding the influence
ED
of the field conditions and properties on 𝑅𝑎𝑆 is shown in Fig. 17. The analysis shows that 𝑘, salinity, 𝜑 and 𝐻 have the highest influence on 𝑅𝑎𝑆 , as indicated by the high absolute values
PT
of their mean EE. However, all the studied parameters display a statistically significant non-
CE
zero impact on 𝑅𝑎𝑆 as they fall below the dashed lines that mark the standard error of the mean, 𝑆𝐷/𝑟 0.5 (see figure caption). Furthermore, high values of SD EE indicate the presence
AC
of non-linearity or interaction effects, and the high influence of salinity on 𝑅𝑎𝑆 , although not visible in the local sensitivity analysis, becomes evident when the average parameter influence over the entire parameter range is considered.
4. Discussion When investigating the prerequisites for enhanced solubility trapping during geological CO2 storage under broad ranges of conditions, we found that the driving force for convective
ACCEPTED MANUSCRIPT mixing ceases under certain thermodynamic and salinity conditions; e.g., under certain conditions, CO2 dissolution into brine causes a density decrease instead of an increase, which concurs with the findings of Lu et al. [74]. Even if the density of the CO2-saturated brine is lighter than the underlying resident brine, this value is still heavier than that of the overlying CO2-rich phase, which leads to a stable density gradient in the formation as long as the CO2-
point, a risk for upward buoyancy-driven flow may exist.
CR IP T
rich phase remains in the formation. However, if lighter CO2-saturated brine reaches a spill
In investigating the relative solutal and direct thermal contributions to the onset of instability
AN US
under broad ranges of conditions, ∆𝜌𝑇 and 𝑅𝑎 𝑇 were found to be highly dependent on the ∇𝑇𝑔𝑒𝑜 and 𝐻 but were generally smaller than ∆𝜌𝑆 and 𝑅𝑎𝑆 , except under certain thermodynamic conditions at high salinity (Fig. 6 and Fig. 7). However, the absolute magnitudes of 𝑅𝑎𝑆 and 𝑅𝑎 𝑇 at these high salinity conditions are both small, resulting in a low
M
𝑅𝑎 that falls slightly above, within or below the 𝑅𝑎𝐶 range. The direct thermal contribution
ED
during these conditions can be quantified by the change in critical depth (depth at which 𝑅𝑎 = 4𝜋 2 ) between when 𝑅𝑎 𝑇 is considered in Eq. 1 and when it is not. Accounting for the direct
PT
thermal contribution leads to an increase in critical depth. For the cold and warm basin conditions (6 m NaCl, 𝜑=0.3, 𝐻=100 m) that are shown in Fig. 9, the critical depths increase
CE
by ~44 m and ~54 m, respectively. Thus, the direct thermal impact, i.e., the increase in the
AC
storage depth for which a prerequisite for density-driven instability exists, is limited to at most a few tens of meters. Several previous studies [8,46,51] have concluded that the thermal contribution to convection is negligible compared to that of the solutal. However, we found that clear, non-negligible, indirect thermal contributions to the onset of instability and convective mixing exist according to the results of the conducted trend and sensitivity analyses, although the direct thermal contribution generally is small. The trends in 𝑅𝑎 (Fig. 9), 𝑅𝑎𝑆 (Fig. 10), 𝑡𝐶 (Fig. 11) and 𝐹 (Fig. 14) were shown to be affected by the basin type
ACCEPTED MANUSCRIPT (i.e., the ∇𝑇𝑔𝑒𝑜 ). The local sensitivity analysis (Fig. 16) identified that the ∇𝑇𝑔𝑒𝑜 influences the 𝑥𝑆 , 𝜇, 𝐷, ∆𝜌𝑆 , 𝑅𝑎𝑆 , 𝐹, 𝑡𝐶 and the convective regime duration. In addition, the global sensitivity analysis confirmed that the ∇𝑇𝑔𝑒𝑜 has a statistically significant non-zero impact on 𝑅𝑎𝑆 .
CR IP T
Investigating the prerequisites for density-driven instability under broad ranges of field conditions (Fig. 8) indicated that combinations of high temperature (≥ 80 °C), moderate-high salinity (≥ 4 m NaCl) and low 𝑘 (≤ 10-14 m2) disfavour density-driven instability. 𝑘 has a
major influence on the magnitude of 𝑅𝑎, which increases with 𝑘. This was also confirmed by
AN US
the local and global sensitivity analyses (Figs. 16 and 17) in this study. However, this study also shows that the trends in 𝑅𝑎 over the broad ranges of thermodynamic and salinity conditions follow the trends in ∆𝜌𝑆 , namely, decreasing with temperature and salinity and increasing with pressure. Thus, deep saline aquifers with relatively low temperature and
M
salinity and high 𝑘 provide the optimal conditions for the onset of density-driven instability,
ED
such as the conditions that were reported for the North Sea reservoir [8]. Indeed, the result from the site-specific 𝑅𝑎𝑆 comparison showed that the Sleipner site has exceptionally good
PT
prerequisites for density-driven instability compared to the other sites that were studied.
CE
Furthermore, cold basin conditions are more favourable to the onset of density-driven instability than warm basin conditions, rendering higher 𝑅𝑎 (Fig. 9) and shorter 𝑡𝐶 (Fig. 11)
AC
for the same depth. Thus, the lowest depth at which the prerequisites for instability exist is shallower for warm basin conditions than for cold basin conditions (Figs. 9 and 10), given the same salinity. Consequently, as this depth decreases with increasing salinity, the choice of storage depth is especially important for CO2 storage in warm basins with high salinity. The general trend that was found in this study for all basin types, salinities, and 𝑘 is that the largest 𝑅𝑎 is obtained at temperature and pressure conditions that correspond to relatively
ACCEPTED MANUSCRIPT shallow depths around 800 m and thereafter decreases with depth. At the same time, the 𝑡𝐶 was found to increase with storage depth for moderate-high salinities (≥3 m NaCl), but this trend ceases with decreasing salinity (Fig. 11), which indicates that CO2 storage at shallower depths favours the onset of density-driven instability and short 𝑡𝐶 , especially within highly saline aquifers. In the study by Bachu [71], the optimal storage depths for free phase CO2
CR IP T
based on capacity and drilling costs were concluded to be 800-1000 m for cold basins and 1500-2000 m for warm basins. For cold basins, the suggested depth range also coincides with the largest 𝑅𝑎. Accordingly, this depth range is optimal from both a capacity (free phase
storage) and onset of instability (enhanced solubility trapping) perspective. For warm basins,
AN US
the optimal depth that was mentioned by Bachu [71] does not coincide with the depth for the largest 𝑅𝑎. The values of 𝑅𝑎 that are found in this depth range only exceed the 𝑅𝑎𝐶 range for salinities <4 m NaCl. These results indicate that convection might not set in at this depth
M
range at higher salinities, instead favouring the onset of instability and short 𝑡𝐶 to store CO2 at shallow depths. However, these shallow depth conditions do not necessarily create the highest
ED
steady convective flux, 𝐹. As shown by the result from the local sensitivity analysis (Fig. 16)
PT
and trend analyses, 𝐹 is highly dependent on the 𝑘 (comparing Figs. 13 and 14), the salinity (Figs. 14 and 16), the amount of exchanges with the capillary transition zone (the upper
CE
boundary permeability, Fig. 13), the basin type (i.e., ∇𝑇𝑔𝑒𝑜 ) and the storage depth (Fig. 14). Both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws show that 𝐹 follows the salinity-
AC
dependent trend of the ∆𝜌𝑆 /𝜇 ratio with depth within an aquifer. However, at moderate to high salinities (≥3 m NaCl), the thermodynamic conditions that benefit both 𝐹 and 𝑡𝐶 coincide at shallow depth if the prerequisites for instability exist. In addition, the results (Figs. 10 and 14) imply that the storage sites with the highest 𝑅𝑎𝑆 do not necessarily have the highest 𝐹.
ACCEPTED MANUSCRIPT The impact of salinity, basin type and storage depth on the size of 𝐹 and the approximate onset and duration of the steady convective regime is compelling (Fig. 15), i.e., the time scale of solubility trapping is substantially affected by the prevailing thermodynamic and salinity conditions of the storage site and is thus of great concern to risk assessment during site characterization.
CR IP T
The findings of this study can be applied in a straight-forward manner to estimate the 𝑅𝑎𝑆 of a specific site or to obtain a greater understanding of which field conditions promote densitydriven instability and convection and enhance solubility trapping. Site-specific estimates of 𝑅𝑎𝑆 can be obtained by using the procedure that was presented in section 3.6 (applying the
AN US
proxy-equations in Appendix A or the result in Fig. 12) without running the developed analysis tool, making it possible to easily evaluate and compare storage locations’
prerequisites for the onset of instability, i.e., enhanced solubility trapping, as part of a site
M
characterization effort. The findings of this study could also be used for future studies together with global field data to map storage formations with conditions that are favourable
ED
for enhanced solubility trapping.
PT
The present study considers density-driven instability and convection in a homogenous isotropic porous layer, not convection that is driven by changes in, e.g., 𝜇 from CO2
CE
dissolution or temperature differences. The direct thermal contribution to convection can be
AC
expected to increase when thermal effects on 𝜇 are also considered. The effect of dissolved CO2 on 𝜇 in aqueous NaCl solutions was previously studied [75,76] and a viscosity model was presented [77]. Neglecting the aqueous CO2 concentration dependency of 𝜇 will result in a slight underestimation of the aqueous phase 𝜇 (increasing with aqueous CO2 concentration), a disregard of an indirect effect on 𝐷 and a slight overestimation of the 𝑅𝑎𝑆 . This was also confirmed by calculations of RaS (results not shown) that were conducted with the abovementioned viscosity model incorporated into the analysis tool. However, including the effect
ACCEPTED MANUSCRIPT of aqueous CO2 concentration on 𝜇 only slightly affects the values of 𝑅𝑎𝑆 , while the trend in 𝑅𝑎𝑆 remains the same, which confirms that the influence of thermodynamic and salinity conditions on 𝜇 are stronger than that of the aqueous CO2 concentration. Furthermore, NaCl was selected to represent the salinity in this study because Na+ is one of
CR IP T
the most commonly cations that are found in groundwater [78]. However, other ions may be present in the formation fluid or dominate at certain sites. These ions influence CO2 solubility; both the ionic strength and the composition (valence and size) of the ions determine the
salting-out effect. Thus, brine that contains divalent cations (e.g., MgCl2 or CaCl2 solutions)
AN US
instead of monovalent cations (e.g., NaCl or KCl solutions) exhibits a lower CO2 solubility at the same concentration [79]. The properties of the aqueous phase, e.g., 𝜇, are also affected by ions. The mentioned effects are expected to influence density-driven convection. Further research could aim to quantify these effects over broad ranges of field conditions and the
M
effect of co-injected impurities on density-driven convection, which can be expected to affect,
PT
5. Conclusions
ED
e.g., the solubility and thereby the induced density changes driving convection.
This study provides an overview of the field conditions for which the prerequisites for
CE
density-driven instability and convection, which enhance solubility trapping, exist. An analysis was conducted over broad thermodynamic, salinity and permeability conditions that
AC
are relevant for CO2 storage in deep saline aquifers. New insight into the trends of 𝑅𝑎, 𝑡𝐶 and 𝐹 over these broad ranges of field conditions is provided and the key parameters that control these are identified. For this purpose, a new analysis tool for comprehensive calculations of 𝑅𝑎, 𝑡𝐶 and 𝐹 was developed to account for the temperature, pressure and salinity dependency of ∆𝜌𝑆 , ∆𝜌𝑇 , 𝜇, 𝐷 and 𝜅 and their combined effect on density-driven instability. These dependences are usually disregarded or simplified but, as this study shows, are of great
ACCEPTED MANUSCRIPT importance to understand the trends of 𝑅𝑎, 𝑡𝐶 and 𝐹 under field conditions. In addition, the relative influence of field parameters, e.g. storage depth and basin type (i.e., geothermal gradient), were analysed through local and global sensitivity analyses. In our study of the conditions that favour and disfavour density-driven instability and convective mixing, we found that the magnitude of 𝑅𝑎 is highly influenced by permeability
CR IP T
but that the trends of 𝑅𝑎 follow the trends of ∆𝜌𝑆 , i.e., decreasing with temperature and salinity and increasing with pressure.
The direct thermal contribution (∆𝜌𝑇 ) to density-driven instability was also found to be small
AN US
compared to the solutal contribution (∆𝜌𝑆 ), except at high salinities. Contrary to previous studies, this study shows that the geothermal gradient has a non-negligible effect on densitydriven instability and convective mixing when accounting for both direct and indirect thermal effects. Cold basin conditions generally render higher 𝑅𝑎 compared to warm. However, the
M
largest 𝑅𝑎 for all basin types and salinities is obtained for thermodynamic conditions that
ED
correspond to relatively shallow depths around 800 m, which then decreases with depth. This indicates that shallow depth CO2 storage, especially at high salinities, favours the onset of
PT
density-driven instability and short 𝑡𝐶 . Shallow depth conditions do not, however, necessarily
CE
render the largest 𝐹, as shown in a trend analysis of both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws for 𝐹; instead, the optimal storage depth for large 𝐹 depends on the salinity,
AC
owing to the ratio ∆𝜌𝑆 /𝜇. Altogether, the salinity, basin type and storage depth highly affect the onset of density-driven instability, 𝑡𝐶 , 𝐹 and the duration of the convective regime, subsequently altering the solubility trapping time scale. In addition, a straight-forward and efficient procedure to estimate site-specific 𝑅𝑎𝑆 was presented. 𝑅𝑎𝑆 may be obtained for any site by means of the dynamic factor (∆𝜌𝑆 /𝜇𝐷), which
ACCEPTED MANUSCRIPT is comprehensively calculated by the developed analysis tool. The estimated 𝑅𝑎𝑆 can serve as an indicator for density-driven instability prerequisites at a site when evaluating its suitability as a storage location.
CR IP T
Acknowledgments The authors would like to thank Shide Mao at the China University of Geosciences for
providing the FORTRAN source code to his viscosity model [59] and Nicholas Spycher at Lawrence Berkeley National Laboratory for providing the FORTRAN source code to his
AN US
solubility model [48,49]. The research that led to these results has received funding from the European Community's 7th Framework Programme MUSTANG and PANACEA projects (grant agreement numbers 227286 and 282900) and Swedish Research Council (VR). The
M
authors would like to thank the anonymous reviewers for their reviews and constructive
ED
comments, which improved this manuscript. Appendix A
PT
Proxy-equations for the dynamic factors are shown in Fig. 12.
CE
Dynamic factor = (𝑎1 + 𝑎2 ∙ T 2 + 𝑎3 ∙ T 3 + 𝑎4 ∙ T 4 + 𝑎5 ∙ 𝑃 + 𝑎6 ∙ P2 + 𝑎7 ∙ P3 + 𝑎8 ∙ P4 + 𝑎9 ∙ T ∙ P + 𝑎10 ∙ 𝑇 2 ∙ 𝑃 + 𝑎11 ∙ 𝑇 ∙ 𝑃2 ) ∙ 1012
Here, the pressure, P, is in bars and the temperature, T, is in ˚C.
AC
The coefficients are given in Table A.
ACCEPTED MANUSCRIPT References [1] IPCC. IPCC Special Report on carbon dioxide capture and storage. Cambridge University Press, Cambridge (2005). [2] Garcia JE. Density of aqueous solutions of CO2. LBNL-49023. Lawrence Berkeley
CR IP T
National Laboratory, Berkeley (2001). [3] Pruess K, Zhang K. Numerical modeling studies of the dissolution-diffusion-convection process during CO2 storage in saline aquifers. LBNL-1243E. Lawrence Berkeley National
AN US
Laboratory, Berkeley (2008).
[4] Ennis-King J, Paterson L. Rate of dissolution due to convective mixing in the underground storage of carbon dioxide. 6th International Conference on Greenhouse Gas Control
M
Technologies 1: 507-510 (2003).
[5] Yang C, Gu Y. Accelerated mass transfer of CO2 in reservoir brine due to density-driven
PT
2430-2436 (2006).
ED
natural convection at high pressures and elevated temperatures. Ind Eng Chem Res 45(8):
[6] Hassanzadeh H, Pooladi-Darvish M, Keith DW. Modelling of convective mixing in CO2
CE
storage. J Can Petrol Technol 44(10): 43-50 (2005).
AC
[7] Weir GJ, White SP, Kissling WM. Reservoir storage and containment of greenhouse gases. Transp Porous Med 23: 37-60 (1996). [8] Lindeberg E, Wessel-Berg D. Vertical convection in an aquifer column under a gas cap of CO2. Energ Convers Manage 38: Suppl., S229-S234 (1997).
ACCEPTED MANUSCRIPT [9] Farajzadeh R, Salimi H, Zitha PLJ, Bruining H. Numerical simulation of density-driven natural convection in porous media with application for CO2 injection projects. Int J Heat Mass Tran 50: 5054-5064 (2007). [10] Kneafsey TJ, Pruess K. Laboratory flow experiments for visualizing carbon dioxide-
CR IP T
induced, density-driven brine convection. Transp Porous Med 82: 123-139 (2010). [11] Tilton N, Riaz A. Nonlinear stability of gravitationally unstable, transient, diffusive boundary layers in porous media. J Fluid Mech 745: 251-278 (2014).
[12] Elenius MT, Johannsen K. On the time scales of nonlinear instability in miscible
AN US
displacement porous media flow. Comput Geosci 16: 901-911 (2012).
[13] Bestehorn M, Firoozabadi, A. Effect of fluctuations on the onset of density-driven
M
convection in porous media. Phys Fluids 24: 114102 (2012).
[14] Riaz A, Hesse M, Tchelepi HA, Orr Jr. FM. Onset of convection in a gravitationally
ED
unstable diffusive boundary layer in porous media. J Fluid Mech 548: 87-111 (2006).
PT
[15] Hassanzadeh H, Pooladi-Darvish M, Keith DW. Stability of a fluid in a horizontal saturated porous layer: effect of non-linear concentration profile, initial, and boundary
CE
conditions. Transp Porous Med 65: 193-211 (2006).
AC
[16] Rapaka S, Chen S, Pawar RJ, Stauffer PH, Zhang D. Non-modal growth of perturbations in density-driven convection in porous media. J Fluid Mech 609: 285-303 (2008). [17] Xu X, Chen S, Zhang D. Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv Water Resour 29: 397-407 (2006). [18] Ennis-King J, Paterson L. Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations. SPE J, September 2005: 349- 356 (2005).
ACCEPTED MANUSCRIPT [19] Ennis-King J, Preston I, Paterson L. Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys Fluids 17: 084107 (2005). [20] Rapaka S, Pawar RJ, Stauffer PH, Zhang D, Chen S. Onset of convection over a transient base-state in anisotropic and layered porous media. J Fluid Mech 641: 227-244 (2009).
driven convection. Transp Porous Med 72: 241-253 (2008).
CR IP T
[21] Hong JS, Kim MC. Effect of anisotropy of porous media on the onset of buoyancy-
[22] Hassanzadeh H, Pooladi-Darvish M, Keith DW. Scaling behavior of convective mixing,
AN US
with application to geological storage of CO2. AIChE J 53(5): 1121-1131 (2007).
[23] Azin R, Raad SMJ, Osfouri S, Fatehi R. Onset of instability in CO2 sequestration into saline aquifer: scaling relationship and the effect of perturbed boundary. Heat Mass Transfer
M
49: 1603-1612 (2013).
[24] Green CP, Ennis-King J. Effect of vertical heterogeneity on long-term migration of CO2
ED
in saline formations. Transp Porous Med 82: 31-47 (2010).
PT
[25] Pau GSH, Bell JB, Pruess K, Almgren AS, Lijewski MJ, Zhang K. High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers. Adv
CE
Water Resour 33(4): 443-455 (2010).
AC
[26] Chen C, Zeng L, Shi L. Continuum-scale convective mixing in geological CO2 sequestration in anisotropic and heterogeneous saline aquifers. Adv Water Resour 53: 175-187 (2013).
[27] Ennis-King J, Paterson L. Coupling of geochemical reactions and convective mixing in the long-term geological storage of carbon dioxide. Int J Greenhouse Gas Control 1: 86-93 (2007).
ACCEPTED MANUSCRIPT [28] Javaheri M, Abedi J, Hassanzadeh H. Onset of convection in CO2 sequestration in deep inclined saline aquifers. J Can Petrol Technol 48(8): 22-27 (2009). [29] Slim AC, Ramakrishnan TS. Onset and cessation of time-dependent, dissolution-driven convection in porous media. Phys Fluids 22: 124103. DOI:10.1063/1.3528009 (2010).
CR IP T
[30] Elenius MT, Nordbotten, JM, Kalisch, H. Effects of a capillary transition zone on the stability of a diffusive boundary layer. IMA J Appl Math 77: 771-787 (2012).
[31] Emami-Meybodi H, Hassanzadeh H. Two-phase convective mixing under a buoyant
AN US
plume of CO2 in deep saline aquifers. Adv Water Resour 76: 55-71 (2015).
[32] Neufeld JA, Hesse MA, Riaz A, Hallworth MA, Tchelepi HA, Huppert HE. Convective dissolution of carbon dioxide in saline aquifers. Geophys Res Lett 37: L22404 (2010).
M
[33] Backhaus S, Turitsyn K, Ecke RE. Convective instability and mass transport of diffusion layers in a Hele-Shaw geometry. Phys Rev Lett 106(10): 104501 (2011).
ED
[34] Zhang W, Li Y, Omambia AN. Reactive transport modeling of effects of convective
PT
mixing on long-term CO2 geological storage in deep saline formations. Int J Greenhouse Gas Control 5: 241-256 (2011).
CE
[35] Hidalgo JJ, MacMinn CW, Juanes R. Dynamics of convective dissolution from a
AC
migrating current of carbon dioxide. Adv Water Resour 62: 511-519 (2013). [36] Szulczewski ML, Hesse MA, Juanes R. Carbon dioxide dissolution in structural and stratigraphic traps. J Fluid Mech 736: 287-315 (2013). [37] Slim AC, Bandi MM, Miller JC, Mahadevan L. Dissolution-driven convection in a HeleShaw cell. Phys Fluids 25: 024101 (2013).
ACCEPTED MANUSCRIPT [38] Slim AC. Solutal-convection regimes in a two-dimesional porous medium. J Fluid Mech 741: 461-491 (2014). [39] Hesse MA. Mathematical modeling and multiscale simulation for CO2 storage in saline aquifers. PhD, Dept. of Energy Resources Engineering, Stanford University (2008).
media. Phys Rev Lett 109: 264503 (2012).
CR IP T
[40] Hidalgo JJ, Fe J, Cueto-Felgueroso L, Juanes R. Scaling of convective mixing in porous
[41] Ghesmat K, Hassanzadeh H, Abedi J. The impact of geochemistry on convective mixing in a gravitationally unstable diffusive boundary layer in porous media: CO2 storage in saline
AN US
aquifers. J Fluid Mech 673: 480-512 (2011).
[42] Islam A, Korrani AKN, Sepehrnoori K, Patzek T. Effects of geochemical reaction on
Mass Tran 77: 519-528 (2014).
M
double diffusive natural convection of CO2 in brine saturated geothermal reservoir. Int J Heat
ED
[43] Farajzadeh R, Ranganathan P, Zitha PLJ, Bruining J. The effect of heterogeneity on the character of density-driven natural convection of CO2 overlying a brine layer. Adv Water
PT
Resour 34: 327-339 (2011).
CE
[44] Ranganathan P, Farajzadeh R, Bruining H, Zitha PLJ. Numerical simulation of natural convection in heterogeneous porous media for CO2 geological storage. Transp Porous Med
AC
95(1): 25-54 (2012).
[45] Kong X-Z, Saar MO. Numerical study of the effects of permeability heterogeneity on density-driven convective mixing during CO2 dissolution storage. Int J Greenhouse Gas Control 19: 160-173 (2013).
ACCEPTED MANUSCRIPT [46] Islam AW, Lashgari HR, Sephernoori K. Double diffusive natural convection of CO2 in a brine saturated geothermal reservoir: study of non-modal growth of perturbations and heterogeneity effects. Geothermics 51: 325-336 (2014). [47] Duan Z, Sun R. An improved model calculating CO2 solubility in pure water and aqueous
(2003).
CR IP T
NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chem Geol 193(3-4): 257-271
[48] Spycher N, Pruess K, Ennis-King J. CO2-H2O mixtures in the geological sequestration of CO2. I. Assessment and calculation of mutual solibilities from 12 to 100˚C and up to 600 bar.
AN US
Geochim Cosmochim Acta 67(16): 3015-3031 (2003).
[49] Spycher N, Pruess K. CO2-H2O mixtures in the geological sequestration of CO2. II. Partitioning in chloride brines at 12-100˚C and up to 600 bar. Geochim Cosmochim Acta
M
69(13): 3309-3320 (2005).
ED
[50] Loodts V, Rongy L, De Wit A. Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions. CHAOS 24: 043120
PT
(2014).
CE
[51] Javaheri M, Abedi L, Hassanzadeh H. Linear stability analysis of double-diffusive convection in porous media, with application to geological storage of CO2. Transp Porous
AC
Med 84: 441-456 (2010). [52] Nield DA, Bejan A. Convection in porous media, 3rd edition. Springer (2006). [53] Kumar A, Ozah R, Noh M, Pope GA, Bryant S, Sepehrnoori K, et al., Reservoir simulation of CO2 storage in deep saline aquifers. SPE J 10(3): 336-348 (2005).
ACCEPTED MANUSCRIPT [54] Li D, Graupner BJ, Bauer S. A method for calculating the liquid density for the CO2H2O-NaCl system under CO2 storage condition. Energy Procedia 4: 3817-3824 (2011). [55] Mao S, Duan Z. The P,V,T,x properties of binary aqueous chloride solutions up to T=573 K and 100 MPa. J Chem Thermodyn 40(7): 1046-1063 (2008).
100 MPa. Energ Fuel 22(3): 1666-1674 (2008).
CR IP T
[56] Duan Z, Hu J, Li D, Mao S. Densities of the CO2-H2O-NaCl systems up to 647 K and
[57] Wagner W, Cooper JR, Dittmann A, Kijima J, Kretzschmar HJ, Kruse A, et al., The IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. J
AN US
Eng Gas Turbines Power 122(1): 150-182 (2000).
[58] Archer DG, Wang P. The dielectric constant of water and Debye-Hückel limiting law
M
slopes. J Phys Chem Ref Data 19(2): 371-411 (1990).
[59] Mao S, Duan Z. The viscosity of aqueous alkali-chloride solutions up to 623 K, 1,000
ED
bar, and high ionic strength. Int J Thermophys 30(5): 1510-1523 (2009).
PT
[60] Adams JJ, Bachu S. Equations of state for basin geofluids: algorithm review and intercomparison for brines. Geofluids 2(4): 257-271 (2002).
CE
[61] Wilke CR, Chang P. Correlation of diffusion coefficients in dilute solutions. AIChE J
AC
1(2): 264-270 (1955).
[62] Cussler EL. Diffusion: mass transfer in fluid systems. 3rd edition, Cambridge University Press, Cambridge (2009). [63] Somerton WH. Thermal properties and temperature-related behavior of rock/fluid systems. Developments in petroleum science, 37. Elsevier, Amsterdam (1992).
ACCEPTED MANUSCRIPT [64] Tikhomirov V. Conductivity of rocks and their relationship with density, saturation and temperature. Neft Khoz 46(4): 36 (1968). [65] Lee Y, Deming D. Evaluation of thermal conductivity temperature corrections applied in terrestrial heat flow studies. J Geophys Res 103(B2): 2447-2454 (1998).
CR IP T
[66] Finsterle S. iTOUGH2 User’s Guide. LBNL-40040. Lawrence Berkeley National Laboratory, Berkeley (2007).
[67] Morris M D. Factorial sampling plans for preliminary computational experiments.
AN US
Technometrics 33(2): 161-174 (1991)
[68] Doherty J, PEST: Model-independent parameter estimation. Watermark numerical computing, Brisbane, Australia (2008), http://pesthomepage.org/
M
[69] Finsterle S, Zhang Y. Solving iTOUGH2 simulation and optimization problems using the PEST protocol. Environmental Modelling & Software 26: 959-968 (2011).
ED
[70] Bachu S. Sequestration of CO2 in geological media in response to climate change: road
PT
map for site selection using the transform of the geological space into the CO2 phase space. Energ Convers Manage 43(1): 87-102 (2002).
CE
[71] Bachu S. Screening and ranking of sedimentary basins for sequestration of CO2 in
AC
geological media in response to climate change. Environ Geol 44(3): 277-289 (2003). [72] Michael K, Allinson G, Golab A, Sharma S, Shulakova V. CO2 storage in saline aquifers II - experience from existing storage operations. Energy Procedia. 1: 1973-1980 (2009). [73] Trémosa J, Castillo C, Quang Vong C, Kervévan C, Lassin A, Audigane P. Long-term assessment of geochemical reactivity of CO2 storage in highly saline aquifers: Application to Ketzin, In Salah and Snøhvit storage sites. Int J Greenhouse Gas Control 20: 2-26 (2014).
ACCEPTED MANUSCRIPT [74] Lu C, Han WS, Lee S-Y, McPherson BJ, Lichtner PC. Effects of density and mutual solubility of a CO2-brine system on CO2 storage in geological formations: "Warm" vs. "cold" formations. Adv Water Resour 32(12): 1685-1702 (2009). [75] Bando S, Takemura F, Nishio M, Hihara E, Akai M. Viscosity of aqueous NaCl solutions
(2004).
CR IP T
with dissolved CO2 at (30 to 60) °C and (10 to 20) MPa. J Chem Eng Data 49: 1328-1332
[76] Fleury M, Deschamps H. Electrical conductivity and viscosity of aqueous NaCl solutions with dissolved CO2. J Chem Eng Data 53: 2505-2509 (2008).
AN US
[77] Islam AW, Carlson ES. Viscosity models and effects of dissolved CO2. Energ Fuel 26: 5330-5336 (2012).
[78] Freeze RA, Cherry JA. Groundwater. Upper Saddle River, Prentice Hall, Inc. (1979). [79] Zerai B, Saylor B Z, Allen DE. Quantification of CO2 trapping and storage capacity in
M
the subsurface: uncertainty due to solubility models. In: McPherson BJ, Sundquist ET,
ED
editors. Carbon sequestration and its role in the global carbon cycle. Geophysical Monograph
AC
CE
PT
Series 183, American Geophysical Union. p. 249-260 (2009).
ACCEPTED MANUSCRIPT
M
AN US
CR IP T
Figure captions
Fig. 1. Schematic of a homogenous, isotropic porous layer. The symbols and abbreviations are
AC
CE
PT
ED
defined in the nomenclature.
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 2. Schematic of the developed analysis tool and its overall structure and submodels.
AN US
CR IP T
ACCEPTED MANUSCRIPT
∆𝜌𝑆 𝐷
+
∆𝜌𝑇 𝑔𝑘𝐻
M
Fig. 3. Comparison of S-numbers, 𝑆 = (
𝜅
)
𝜇
, from the study by Lindeberg and
Wessel-Berg [8] and those that were calculated by the developed analysis tool, which account
AC
CE
PT
ED
for constant and varying ∆𝜌𝑆 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 4. Solubility of CO2, 𝑥𝑆 , at different temperatures and pressures for salinities of (a) 1 m
M
NaCl, (b) 4 m NaCl and (c) 6 m NaCl. The solid and dashed lines illustrate the results of the
AC
CE
PT
ED
solubility models by Duan and Sun [47] and Spycher and Pruess [49], respectively.
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 5. Density changes from CO2 dissolution, ∆𝜌𝑆 %, at different temperatures and pressures
M
for salinities of (a) 1 m NaCl, (b) 4 m NaCl and (c) 6 m NaCl. The solid and dashed lines
AC
CE
PT
[49], respectively.
ED
illustrate the results of the solubility models by Duan and Sun [47] and Spycher and Pruess
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 6. Buoyancy ratios, 𝑁 = ∆𝜌𝑆 /∆𝜌𝑇 , at various depths for different geothermal gradients at
AC
CE
PT
ED
M
salinities of (a) 0 m NaCl and (b) 6 m NaCl. Hydrostatic pressure is assumed.
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 7. The ratio 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 at various depths for different geothermal gradients at salinities
AC
CE
PT
ED
M
of (a) 0 m NaCl and (b) 6 m NaCl. Hydrostatic pressure is assumed.
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
Fig. 8. 𝑅𝑎 for different temperature, pressure, salinity and permeability conditions. The left column shows the results for k=10-12 m2, the middle column for k=10-14 m2 and the right column for k=10-15 m2. The coloured curves indicate different pressures: black - 7.85 MPa, green - 30 MPa, and blue - 60 MPa. The solid and dashed lines illustrate the results of the solubility models by Duan and Sun [48] and Spycher and Pruess [49], respectively. The solid red lines indicate the 𝑅𝑎𝐶 range: 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 9. 𝑅𝑎 for cold basin (black lines) and warm basin (magenta lines) conditions with ∇𝑇𝑔𝑒𝑜 = 25 °C/km and 50 °C/km, respectively. 𝑇𝑠𝑢𝑟𝑓 =10 °C and hydrostatic pressure is assumed. 𝑅𝑎
ED
values are shown for different permeabilities and salinities of 0 m NaCl, 1 m NaCl, 4 m NaCl
AC
CE
PT
and 6 m NaCl. The solid red lines indicate the 𝑅𝑎𝐶 range: 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 10. 𝑅𝑎𝑆 at different salinities and various depths. To the left: for a cold basin with ∇𝑇𝑔𝑒𝑜 = 25 °C/km. To the right: for a warm basin with ∇𝑇𝑔𝑒𝑜 = 50 °C/km. For both cases, 𝑇𝑠𝑢𝑟𝑓 =10 °C
ED
and k=10-13 m2, and hydrostatic pressure is assumed. The solid red lines indicate the 𝑅𝑎𝐶
AC
CE
PT
range: 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 11. The 𝑡𝐶 at different salinities and various depths. To the left: for a cold basin with
AC
CE
PT
ED
𝑇𝑠𝑢𝑟𝑓 =10 °C and k=10-13 m2.
M
∇𝑇𝑔𝑒𝑜 =25 °C/km. To the right: for a warm basin with ∇𝑇𝑔𝑒𝑜 =50 °C/km. For both cases,
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
CE
Fig. 12. Dynamic factor, ∆𝜌𝑆 /𝜇𝐷𝐶𝑂2 , for different temperatures, pressures and salinities. The
AC
arrow indicates the direction of increasing pressure. Pressure increases with 10 MPa for each curve. The first curve corresponds to 10 MPa, the second curve corresponds to 20 MPa and so on up to 60 MPa.
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 13. Comparison of the scaling laws of the steady convective flux. To the left is the non-
M
dimensional mass transfer (𝐹 ∗ ) at zero salinity for a cold basin with ∇𝑇𝑔𝑒𝑜 =25 °C/km,
ED
𝑇𝑠𝑢𝑟𝑓 =10 °C, k=10-12 m2 and a hydrostatic pressure gradient. The dotted lines represent
AC
CE
PT
experimental uncertainty bounds. To the right is the corresponding dimensional flux (𝐹).
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 14. To the left is the dimensional steady convective flux, 𝐹, at three different salinities and various depths within a cold basin (black symbols) and a warm basin (magenta symbols)
ED
with ∇𝑇𝑔𝑒𝑜 =25 °C/km and 50 °C/km, respectively. 𝑇𝑠𝑢𝑟𝑓 =10 °C and k=10-13 m2, and
AC
CE
PT
hydrostatic pressure is assumed. To the right is the corresponding ∆𝜌𝑆 /𝜇.
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 15. Convective regime duration at two different depths and salinities for a cold basin (solid lines) and a warm basin (dashed lines) with ∇𝑇𝑔𝑒𝑜 =25 °C/km and 50 °C/km,
AC
CE
PT
ED
respectively, and a hydrostatic pressure gradient. 𝑇𝑠𝑢𝑟𝑓 =10 °C and k=10-13 m2.
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 16. Scaled sensitivity coefficients, which show the relative influence of different field parameters on 𝑥𝑆 , ∆𝜌𝑆 , 𝐷, 𝜇, 𝑅𝑎𝑆 , 𝐹, 𝑡𝐶 and the convective regime duration. 𝐹 a is based on the
AC
CE
PT
Neufeld et al. [32].
ED
𝑅𝑎𝑆 -independent scaling law by Slim [38]. 𝐹 b is based on the 𝑅𝑎𝑆 -dependent scaling law by
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 17. Mean EE, which shows the relative influence of different field parameters on 𝑅𝑎𝑆 .
M
SD EE conveys non-linearity or interaction effects. The dashed lines represent the mean EE
AC
CE
PT
ED
±2 (SD/𝑟 0.5).
ACCEPTED MANUSCRIPT
Table captions Table 1. 𝑅𝑎𝐶 for a homogenous and isotropic porous medium under different upper boundary conditions. For all cases, the lower boundary is assumed to be impermeable to fluid and heat/solute flow, with a prescribed temperature and/or solute concentration. 𝑅𝑎𝐶 4𝜋 27.10
impermeable to fluid and heat/solute flow, constant 𝑇 and 𝐶𝑆 .a
26.09
partially permeable (0.1) to fluid flowb, constant 𝐶𝑆 .c
23.45
partially permeable (0.5) to fluid flowb, constant 𝐶𝑆 .c
21.76
fully permeable to fluid flowb, constant 𝐶𝑆 .c
impermeable to fluid flow, constant 𝑃, 𝑇 and 𝐶𝑆 .a
CR IP T
Upper boundary conditions
2
AN US
13.8 impermeable to fluid flow, constant 𝑃 and 𝐶𝑆 .d a Based on linear stability analysis for thermally and solutally induced instability with a steady base state (initial temperature and concentration distributions that vary linearly with depth) [52]. b
Continuity of flux and pressure. The relative permeability of the capillary transition zone is given within brackets. c
Numerically derived for solutally (CO2) induced instability with an evolving base state [29].
d
M
Minimum value that is obtained from (energy) stability boundary mapping for solutally (CO2) induced instability with an evolving base state [29,37].
F*
Valid 𝑅𝑎𝑆 range
Num.
impermeable
-
Num.
[25]
impermeable
-
Num.[38]
impermeable
102-5∙104
Num.[30]
impermeable
Method
0.017
CE
0.017 to 0.018 0.017 a 0.02
AC
0.025 0.044
Upper boundary
[39]
PT
convective flux, 𝐹 ∗ .
ED
Table 2. Numerically and experimentally derived scaling laws of the non-dimensional steady
0.075 (0.12 ± 0.03)𝑅𝑎𝑆 −0.16±0.02
Num.
[38]
10 -5∙104
partially permeable (0.6)c
102-5∙104
partially permeable (0.2)
Num.[38] Num.
[30]
Analogue fluid exp.
c
fully permeable [32]
-
2
b
3
10 -106
(0.045 ± 0.025)𝑅𝑎𝑆 −0.24±0.06 Analogue fluid exp.[33] 6∙103-9∙104 a ∗ Duration of the convective regime is approximated by the parameterization 𝐹 = 0.017, 1100 < 𝑡 < 16 ∙ 𝑅𝑎𝑆 by Slim [38], where 𝑡 is the non-dimensional time that is obtained by multiplying the dimensional time with (𝑔𝑘∆𝜌𝑆 )2 (𝜇𝜑)2 𝐷
.
b
Full exchange with the capillary transition zone.
c
Relative permeability of the capillary transition zone.
ACCEPTED MANUSCRIPT
Table 3. Aquifer properties and conditions that were used to calculate 𝑁, 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 , 𝑅𝑎 and 𝐹 for cold and warm basins. 𝑇𝑠𝑢𝑟𝑓 =10 °C, 𝜑=0.3 and 𝐻=100 m.
𝑁 and 𝜑Ra S /RaT
1.0∙10
a
𝛻𝑇𝑔𝑒𝑜 [°C/km]
Salinity [m NaCl]
-
25, 30, 40, 50
0-6 0-6
-15
1.0∙10 -1.0∙10
-12
7.85-60
30-130
50
-15
-12
a
-
25, 50
1.0∙10 -1.0∙10
𝐹-scaling laws 𝐹 trends
T [°C]
a
1.0∙10-12
a
-13
a
1.0∙10 A hydrostatic pressure gradient is assumed.
0-6
25
0
25, 50
0, 1, 3
AN US
Ra Basin types
P [MPa]
-12
CR IP T
𝑘 [m2]
Cases
Table 4. Properties of the CO2 pilot test and storage sites that were used as input parameters to calculate the site-specific 𝑅𝑎𝑆 and 𝑡𝐶 .
𝑘 [m ] 2
Ketzin73
24
80
1.5∙10
𝜑 [-]
0.3
T [°C]
56
P [MPa]
-12
15.2
a
Sleipner73
M
H [m]
Frio73
7.5∙10
-13
0.23
In Salah74
Snøhvit74
20
75
250 5.0∙10
-12
1.0∙10
-14
8.83∙10-13
0.37
0.15
0.15
34
37
96
98
7.3
10.3
18.0
28.5
ED
Parameters
CE
PT
Salinity [m NaCl] 1.598 4.288 0.587 2.712 2.566 a Salinities converted to m NaCl from TDS, taking into account the site-specific thermodynamic conditions.
AC
Table 5. Parameters that were used in the sensitivity analyses. 𝑘 [m2]
Depth [m]
𝜑 [-]
1.0∙10-13
1000
0.3
𝐻 [m]
∇𝑇𝑔𝑒𝑜 [°C/km]
Salinity [m NaCl]
25
1
a
Local analysis initial value
100
Global analysisb range 1.0∙10-15-1.0∙10-11 800-1800 0.1-0.4 10-100 25-50 0-6 a The initial parameter values were perturbed by 10 % and centred finite differences were used. The sensitivity coefficients were scaled with the ratio of the parameter and variable variations, which were assumed to be 10 % of their initial values. b
Parameter ranges were partitioned into 11 equal intervals and 800 paths were used (= 5600 simulations).
ACCEPTED MANUSCRIPT Table 6. Calculated ∆𝜌𝑆 and 𝑅𝑎𝑆 for specific CO2 pilot test and storage sites.
∆𝜌𝑆 [kg/m3] ∆𝜌𝑆 % [%] 𝑅𝑎𝑆 [-] 𝑡𝐶 [days] a
Frio
Ketzin
Sleipner
In Salah
Snøhvit
5.3
2.7
10.1
1.5
2.0
0.51
0.24
0.99
0.14
0.19
3113.9
3748.0
177,702.5
8.6
3870.9
0.3
c
6.3
102.0
23.8
c 𝑡𝐶 b [days] 29.5 480.5 1.3 a Corresponding to permeable upper boundary condition, 𝑐0 = 31, [30].
CR IP T
112.2
b
Corresponding to impermeable upper boundary condition, 𝑐0 = 146, [14].
c
𝑡𝐶 is not calculated since stable conditions prevail, i.e., 𝑅𝑎𝑆 < 𝑅𝑎𝐶 .
AN US
Table A. Coefficients for the proxy-equations. 0.0 7.948136110915495 -0.003441435234038 0.000040935045132 -0.000000139741972 0.018651779590719 -0.000047204955860 0.000000078802468 -0.000000000046039 0.000006178398246 -0.000000060718778 -0.000000034826392 0.99820
0.5 6.707809604348991 -0.002930182594837 0.000034867612722 -0.000000118953377 0.015363104245509 -0.000037208693130 0.000000062337161 -0.000000000036472 0.000000757544439 -0.000000067646659 -0.000000026725659 0.99832
1.0 5.670286672377516 -0.002501054900298 0.000029768688332 -0.000000101467421 0.012661286761162 -0.000029034768917 0.000000048870428 -0.000000000028640 -0.000003567264417 -0.000000071268675 -0.000000020138373 0.99845
m NaCl 𝑎1 𝑎2 𝑎3 𝑎4 𝑎5 𝑎6 𝑎7 𝑎8 𝑎9 𝑎10 𝑎11 𝑅2
1.5 4.801002897984335 -0.002140719463234 0.000025483216074 -0.000000086761195 0.010429841104933 -0.000022323006341 0.000000037804293 -0.000000000022195 -0.000006990943279 -0.000000073206522 -0.000000014734211 0.99858
2.0 4.071765097530080 -0.001837906559726 0.000021879189254 -0.000000074386526 0.008580958668878 -0.000016808787401 0.000000028705856 -0.000000000016889 -0.000009651589099 -0.000000074281849 -0.000000010321768 0.99871
2.5 3.459177590320337 -0.001583147388553 0.000018845297093 -0.000000063964377 0.007044357126234 -0.000012279544781 0.000000021228394 -0.000000000012520 -0.000011663239757 -0.000000075013892 -0.000000006753049 0.99883
3.0
3.5
4.0
ED
PT
CE
AC m NaCl
M
m NaCl 𝑎1 𝑎2 𝑎3 𝑎4 𝑎5 𝑎6 𝑎7 𝑎8 𝑎9 𝑎10 𝑎11 𝑅2
ACCEPTED MANUSCRIPT
2.509486329629760 -0.001187411168486 0.000014129731257 -0.000000047755012 0.004691213083351 -0.000005513538486 0.000000010056460 -0.000000000005972 -0.000014106958036 -0.000000076689295 -0.000000001686914 0.99907
2.142827424385575 -0.001034325676134 0.000012304981938 -0.000000041478565 0.003790733063294 -0.000003017413526 0.000000005938987 -0.000000000003547 -0.000014685146127 -0.000000078002192 -0.000000000009779 0.99917
CR IP T
2.943800700123073 -0.001368518645321 0.000016288158757 -0.000000055176110 0.005763182036722 -0.000008561884080 0.000000015089251 -0.000000000008926 -0.000013121673473 -0.000000075742015 -0.000000003907715 0.99895
AC
CE
PT
ED
M
AN US
𝑎1 𝑎2 𝑎3 𝑎4 𝑎5 𝑎6 𝑎7 𝑎8 𝑎9 𝑎10 𝑎11 𝑅2
Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions Maria Rasmusson , Fritjof Fagerlund , Yvonne Tsang , Kristina Rasmusson , Auli Niemi PII: DOI: Reference:
S0309-1708(15)00195-5 10.1016/j.advwatres.2015.08.009 ADWR 2446
To appear in:
Advances in Water Resources
Received date: Revised date: Accepted date:
27 October 2014 8 August 2015 11 August 2015
Please cite this article as: Maria Rasmusson , Fritjof Fagerlund , Yvonne Tsang , Kristina Rasmusson , Auli Niemi , Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions, Advances in Water Resources (2015), doi: 10.1016/j.advwatres.2015.08.009
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights We develop an instability onset, onset time and steady convective flux-analysis tool. A systematic analysis of prerequisites for enhanced solubility trapping is presented. Indirect thermal effects on enhanced solubility trapping are non-negligible. Shallow storage depths favor density-driven instability onset and short onset times. Salinity determines the storage depth at which the largest convective fluxes occur.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Title: Prerequisites for density-driven instabilities and convective mixing under broad geological CO2 storage conditions
CR IP T
Authors: Maria Rasmussona, Fritjof Fagerlunda, Yvonne Tsangb, Kristina Rasmussona, Auli
Corresponding author: M. Rasmussona
a
M
Affiliations:
AN US
Niemia
Uppsala University – Department of Earth Sciences, Villavagen 16, SE-75236 Uppsala,
ED
Sweden.
PT
+46184717162, [email protected]
CE
+46184717166, [email protected]
AC
+46184717162, [email protected] +46184712263, [email protected] fax number +4618551124 (Department of Earth Sciences) b
Lawrence Berkeley National Laboratory, University of California, Earth Sciences Division, 1
Cyclotron Road, Berkeley, CA 94720-8126, USA. +15104867047, [email protected]
ACCEPTED MANUSCRIPT Abstract Direct atmospheric greenhouse gas emissions can be greatly reduced by CO2 sequestration in deep saline aquifers. One of the most secure and important mechanisms of CO2 trapping over large time scales is solubility trapping. In addition, the CO2 dissolution rate is greatly enhanced if density-driven convective mixing occurs. We present a systematic analysis of the
CR IP T
prerequisites for density-driven instability and convective mixing over the broad temperature, pressure, salinity and permeability conditions that are found in geological CO2 storage. The onset of instability (Rayleigh-Darcy number, 𝑅𝑎), the onset time of instability and the steady
AN US
convective flux are comprehensively calculated using a newly developed analysis tool that accounts for the thermodynamic and salinity dependence on solutally and thermally induced density change, viscosity, molecular and thermal diffusivity. Additionally, the relative influences of field characteristics are analysed through local and global sensitivity analyses.
M
The results help to elucidate the trends of the 𝑅𝑎, onset time of instability and steady
ED
convective flux under field conditions. The impacts of storage depth and basin type (geothermal gradient) are also explored and the conditions that favour or hinder enhanced
PT
solubility trapping are identified. Contrary to previous studies, we conclude that the geothermal gradient has a non-negligible effect on density-driven instability and convective
CE
mixing when considering both direct and indirect thermal effects because cold basin
AC
conditions, for instance, render higher 𝑅𝑎 compared to warm basin conditions. We also show that the largest 𝑅𝑎 is obtained for conditions that correspond to relatively shallow depths, measuring approximately 800 m, indicating that CO2 storage at such depths favours the onset of density-driven instability and reduces onset times. However, shallow depths do not necessarily provide conditions that generate the largest steady convective fluxes; the salinity determines the storage depth at which the largest steady convective fluxes occur.
ACCEPTED MANUSCRIPT Furthermore, we present a straight-forward and efficient procedure to estimate site-specific solutal 𝑅𝑎 that accounts for thermodynamic and salinity dependence.
Concentration [𝑘𝑔𝐶𝑂2 ⁄𝑚3 ]
𝑐
Specific heat capacity [ 𝐽 ⁄(𝑘𝑔 ∙ 𝐾)]
𝑐0
Numerically derived constant [-]
𝐷
Molecular diffusivity [𝑚2 ⁄𝑠]
EE
Mean elementary effect
𝐹
Steady convective flux [𝑘𝑔𝐶𝑂2 ⁄(𝑚2 ∙ 𝑠)]
𝐹∗
Steady convective flux [-]
𝑔
Magnitude of gravitational acceleration [𝑚⁄𝑠 2 ]
𝐻
Layer thickness [𝑚]
𝑘
Permeability [𝑚2 ]
M
ED
PT
CE
Molecular weight [𝑔⁄𝑚𝑜𝑙 ]
AC
̃ 𝑀
AN US
𝐶
CR IP T
Nomenclature
𝑁
Buoyancy ratio [-]
P
Pressure [𝑀𝑃𝑎]
r
Number of paths in Morris OAT method
𝑅𝑎
Rayleigh-Darcy number [-]
ACCEPTED MANUSCRIPT Temperature [˚C (if not otherwise specified)]
𝑡𝑐
Onset time of instability [s]
𝑆
S-number [-]
Salinity
Salinity [𝑚𝑜𝑙𝑁𝑎𝐶𝑙 ⁄𝑘𝑔𝑤 ]
SD
Standard deviation
𝑉̅
Molar volume [𝑐𝑚3 ⁄𝑚𝑜𝑙 ]
𝑥
Solubility [𝑚𝑜𝑙𝐶𝑂2 ⁄𝑘𝑔𝑤 ]
AN US
CR IP T
T
Greek letters
Induced density change [𝑘𝑔⁄𝑚3 ] or [%]
𝛻𝑇𝑔𝑒𝑜
Geothermal gradient [°𝐶 ⁄𝑘𝑚]
𝜃
Empirical constant
𝜅
Thermal diffusivity [𝑚2 ⁄𝑠]
ED
PT
CE
Thermal conductivity of medium [𝑊/(𝑚 ∙ 𝐾)]
AC
𝜆
M
∆𝜌
𝜇
Viscosity [𝑘𝑔⁄(𝑚 ∙ 𝑠)]
𝜌
Density [𝑘𝑔⁄𝑚3 ]
𝜑
Porosity [-]
ACCEPTED MANUSCRIPT
Critical
𝑓
Formation fluid
𝑆
Solutal
𝑠𝑎𝑡
Saturation
𝑠𝑢𝑟𝑓
Surface (temperature)
𝑇
Thermal
𝑤
Water
%
Percentage
1. Introduction
ED
M
AN US
𝑐
CR IP T
Subscripts
PT
CO2 storage in sedimentary basins that contain deep saline aquifers reduces the direct release of greenhouse gases into the atmosphere. The method entails injecting CO2 at ~1-2 km depth
CE
into formations of sufficient porosity, permeability and injectivity. Depending on the
AC
thermodynamic conditions, the injected CO2-rich phase is in a gaseous, liquid or supercritical state. The injected CO2 is retained in the formation by several mechanisms: physical (structural, stratigraphic and hydrodynamic), residual, solubility and mineral trapping [1]. Over large time scales, solubility trapping is one of the most important trapping mechanisms, providing more secure storage compared to free phase CO2 containment [1]. However, CO2 solubility in the aqueous phase (brine) depends on the temperature, pressure and salinity, and the same applies to the accompanying density change in the aqueous phase. Density increases
ACCEPTED MANUSCRIPT of up to 2-3 % are reported following CO2 dissolution into pure water at 10 MPa and up to 110 °C [2], while more moderate density increases up to ~1 % occur in brine at 12-22 MPa and 45-75 °C [3]. These density changes may produce unstable density gradients within aquifers and trigger convective flow, which substantially enhances the CO2 dissolution rate beyond that of diffusion [4,5], thus reducing the ultimate solubility trapping time scale [6].
CR IP T
The connection of density-driven convection to CO2 storage was identified early in numerical [7] and theoretical [8] studies. Enhanced CO2 dissolution rates from density-driven
convection were later verified both experimentally [5,9,10] and numerically [3,4,9]. However, changes in the aqueous phase’s density depend not only on CO2 dissolution but also on the
AN US
geothermal gradient across the formation layer, i.e., density changes are both solutally and thermally induced. Preconditions for density-driven instability can be identified by comparing the Rayleigh-Darcy number, 𝑅𝑎, (Eq. 1, section 2.2), against a critical Rayleigh-Darcy
M
number, 𝑅𝑎𝐶 .
ED
Furthermore, the 𝑅𝑎 is also related to the onset time of instability, 𝑡𝑐 , (Eq. 2, section 2.3), which determines the CO2 solubility trapping time scale. The nonlinear onset of convection
PT
occurs much later than the onset of instability [11] and depends on the initial perturbation [12,13]. The 𝑡𝑐 was studied by using linear, global or non-modal stability analysis for
CE
homogenous [14,15,16] and anisotropic [4,17-21] systems, while the onset of convection was
AC
studied numerically for homogeneous [3,11,22,23] and heterogeneous [24-26] systems. Furthermore, geochemical reactions [27], aquifer inclination [28], and exchanges with the overlying capillary transition zone [29-31] were identified to impact onset time. Several dynamic regimes of the CO2 dissolution process were identified [12,14,22,32-38]. For high 𝑅𝑎, a steady convective flux regime exists [3,25,32,33,38-40] that occurs after the onset of instability and before layer saturation. The steady convective flux, 𝐹, which governs the enhanced CO2 dissolution rate, was quantified both experimentally [32,33,37] and
ACCEPTED MANUSCRIPT numerically [25,30,38,39]. However, discrepancies exist between the obtained scaling laws of 𝐹, particularly with regards to the dependence on 𝑅𝑎. While some studies indicate the nonlinear scaling of 𝐹 with 𝑅𝑎 [32,33], other studies indicate 𝐹 to be 𝑅𝑎-independent [25,3840] but followed by a 𝑅𝑎-dependent shut-down regime [38], thus making the duration of the convective regime 𝑅𝑎-dependent. In addition, 𝐹 was found to increase when exchanges with
CR IP T
the overlying capillary transition zone are considered [30,31]. Other factors that were found to influence density-driven convection are geochemical reactions [34,41,42] and geological heterogeneity [26,43-46]. Heterogeneity affects the flow pattern, changing the flow regime and transport mechanism from convective fingering to dispersive flux with increasing
AN US
heterogeneity [26,43,44]. Thus, the onset of instability is a necessary but insufficient
condition for density-driven convection to occur. However, identifying the field conditions that are the prerequisites for density-driven instability is still paramount.
M
The thermodynamic and salinity conditions control CO2 solubility [47-49] and the
ED
subsequently induced density change, which together with the formation properties, e.g., layer thickness and permeability, determine if the prerequisites for the onset of instability exist. In
PT
addition, increases in salinity were found to delay the onset of instability [34,50], while increases in pressure and decreases in temperature cause destabilization [50]. In earlier
CE
studies, prerequisites for density-driven instability were identified for several sites
AC
[6,8,15,22,51], but the prerequisites for the onset of density-driven instability during CO2 storage have not yet been fully and systematically mapped over the broad range of conditions that are encountered in the field. This study focuses on such systematic trend analyses regarding how 𝑅𝑎, 𝑡𝑐 and 𝐹 evolve over the temperature, pressure, salinity and permeability ranges that are encountered in deep saline aquifers. Additionally, the effects of thermodynamic and salinity conditions on the convective regime duration are investigated.
ACCEPTED MANUSCRIPT Furthermore, the thermal contribution to convection, which is induced by the geothermal gradient, was previously found to be negligible compared to the solutal contribution [8,46,51]. However, we extend the analysis in this study to include both direct and indirect geothermal effects while analysing the influences of the geothermal gradient (i.e., basin type) on the trends in the onset of instability, 𝑡𝑐 and 𝐹.
CR IP T
In addition, the relative influence of different field conditions and properties - such as storage depth, basin type, porosity, permeability, layer thickness and salinity - on density-driven
instability and convective mixing is identified through local and global sensitivity analyses.
AN US
The theoretical findings are discussed in relation to site-specific cases in the form of densitydriven instability prerequisites for CO2 pilot test and storage sites.
For this study, a new analysis tool was developed with the aim to avoid using simplifying
M
assumptions that are commonly made during the estimation of 𝑅𝑎. The tool was used to rigorously calculate 𝑅𝑎, 𝑡𝑐 and 𝐹 by fully accounting for the temperature, pressure and
ED
salinity dependency of CO2 solubility, induced density changes, diffusivities and viscosity.
PT
This approach was possible through the implementation of state-of-the-art solubility, density and viscosity models and algorithms for calculating diffusivities.
CE
These findings give new insight into the field conditions that favour or hinder enhanced
AC
solubility trapping and identify the key field parameters that affect the process; thus, they are expected to be useful during the selection of appropriate CO2 storage sites and for capacity estimation and risk assessment. In addition, a simple and efficient procedure to obtain an accurate estimate of site-specific solutal 𝑅𝑎 is presented.
2. Methods 2.1 Description of the system
ACCEPTED MANUSCRIPT The conceptual model consists of an isotropic, homogenous porous medium of thickness 𝐻 (Fig. 1). Injected CO2 has migrated to the upper part of the system because of buoyancy. Two regions are distinguished: an upper multi-phase region that contains CO2-rich phase and CO2saturated aqueous phase, and an underlying single-phase region that contains aqueous phase with dissolved CO2. The CO2 that is dissolved in the single-phase region induces density
CR IP T
changes in the aqueous phase and, potentially, instability. Several upper boundary conditions are considered: impermeable, partially or fully permeable to flow. These correspond to
different amounts of exchange with the capillary transition zone. In the single-phase region, the CO2 concentration is initially zero (although in some studies assumed to vary linearly with
AN US
depth, see Table 1), while the temperature profile is linear and set by the geothermal gradient.
2.2 Onset of Instability
M
Prerequisites for density-driven instability exist when 𝑅𝑎 > 𝑅𝑎𝐶 . 𝑅𝑎, which is given by Eq. 1,
number, 𝑅𝑎𝑆 [51,52]:
𝑤ℎ𝑒𝑟𝑒
PT
𝑅𝑎 = 𝑅𝑎 𝑇 + 𝑅𝑎𝑆
ED
is the sum of the thermal Rayleigh-Darcy number, 𝑅𝑎 𝑇 , and the solutal Rayleigh-Darcy
𝑅𝑎 𝑇 =
𝑔𝑘∆𝜌𝑇 𝐻 , 𝜇𝜅
𝑅𝑎𝑆 =
𝑔𝑘∆𝜌𝑆 𝐻 𝜇𝜑𝐷
(1)
CE
Here, 𝑔 is the magnitude of gravitational acceleration, 𝑘 is the permeability, 𝐻 is the layer thickness, 𝜇 is the formation fluid viscosity, 𝜑 is the porosity, and ∆𝜌𝑆 and ∆𝜌𝑇 are,
AC
respectively, the solutally and thermally induced density changes. 𝐷 and 𝜅 represent the molecular and thermal diffusivity, respectively. Note that 𝜇, ∆𝜌𝑆 , ∆𝜌𝑇 , 𝐷 and 𝜅 are all temperature, pressure and salinity dependent. Values of 𝑅𝑎𝐶 for different boundary conditions of the system in Fig. 1 are given in Table 1. 𝑅𝑎𝐶 decreases with increasing permeability of the upper boundary. In this study, the 𝑅𝑎𝐶 range 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 was adopted. This 𝑅𝑎𝐶
ACCEPTED MANUSCRIPT range, although 𝑅𝑎𝐶 nearly triples, can be considered narrow when compared to the 𝑅𝑎 range that was studied.
The linear onset time of instability, 𝑡𝑐 , is given by Eq. 2 [4]:
𝑡𝑐 = 𝑐0
1 𝐻2 (𝜇𝜑)2 𝐷 ∙ = 𝑐 0 (𝑅𝑎𝑆 )2 𝐷 (𝑔𝑘∆𝜌𝑆 )2
CR IP T
2.3 Onset time of Instability
(2)
AN US
The value of the numerically derived constant 𝑐0 was reported [3,12,14,15,17-
19,22,23,25,29,30] to range from 𝑐0 =31 [30] to 𝑐0 =5000 [25], depending on the method, initial conditions and stability criteria that were used [22]. In this study, the lower and upper limits of 𝑐0 =31 and 𝑐0 =146, which correspond to full exchange or no exchange with the
ED
2.4 Steady convective flux
M
capillary transition zone [30], were used.
PT
Scaling laws for the non-dimensional steady convective flux, 𝐹 ∗ , are given in Table 2. Both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws are included. Values of 𝐹 ∗ for several
CE
upper boundary conditions with different permeabilities are given. The dimensional steady
AC
convective flux, 𝐹, is given by Eq. 3 [37]: 𝐹 = 𝐹∗ ∙
𝑔𝑘∆𝜌𝑆 ∙ 𝐶𝑠𝑎𝑡 𝜇
(3)
Here, 𝐶𝑠𝑎𝑡 is the mass of solute per unit volume of saturated solution. The onset time of convection depends on the initial perturbation size [12]. A numerically based parameterization that approximates the steady convective regime and its duration was recently presented by Slim [38] and is included in Table 2.
ACCEPTED MANUSCRIPT
2.5 Developed analysis tool A comprehensive 𝑅𝑎-𝑡𝐶 -𝐹-calculation model was developed and coded in MATLAB. See Fig. 2 for the model structure with submodels and work flow. To ensure the accurate and continuous calculation of 𝑅𝑎, 𝑡𝐶 , and 𝐹, the analysis tool incorporates submodels that account
CR IP T
for the temperature, pressure and salinity dependence of CO2 solubility into brine, thermally and solutally induced density changes, the viscosity of the brine and the thermal and
molecular diffusivity. The analysis tool is valid for a broad range of temperatures, 30-130 °C; pressures, 7.85-60 MPa; and salinities, 0-6 m NaCl (mol NaCl/kg H2O). This salinity range,
AN US
which is equivalent to 0-350,000 ppm, was previously applied by Kumar et al. [53], who studied geological CO2 storage trapping mechanisms. 2.5.1 Solubility submodels
M
Two solubility models, Duan and Sun [47] and Spycher and Pruess [48,49], were incorporated
ED
into the analysis tool. Both models are widely used to calculate CO2 solubility under CO2 storage conditions. One of the solubility models, [49], is applicable for temperatures of 12-
PT
100 °C, pressures of 0-60 MPa and salinities of 0-6 m NaCl, while the other model, [47], is valid for temperatures of 0-260 °C, pressures of 0-200 MPa and salinities of 0-4.3 m NaCl.
CE
The use of both solubility models provides a comparison opportunity. The CO2 solubility, 𝑥𝑆 ,
AC
increases with pressure and decreases with ionic strength, while the temperature dependence is non-monotonic. At low pressures, such as <10 MPa, solubility decreases with temperature. At higher pressures, solubility decreases with temperature until a pressure dependent minimum and then increases with temperature [47]. 2.5.2 Density submodel
ACCEPTED MANUSCRIPT A slightly modified version of the aqueous phase density model by Li et al. [54], which accounts for dissolved CO2, was incorporated into the analysis tool. The model by Li et al. [54] combines the density model by Mao and Duan [55] for aqueous NaCl solutions with a similar approach to that of Duan et al. [56] for the density of CO2-H2O-NaCl solutions. Our modifications consisted of adopting the thermodynamic formulation IAPWS97 [57] for pure
CR IP T
water density and calculating the volumetric Debye-Hückel limiting law slope by the linearization of an extensive dataset [58] of tabulated values. Furthermore, the analysis tool calculates ∆𝜌𝑆 as the density difference between the CO2-H2O-NaCl solution and the resident brine. The ∆𝜌𝑇 is calculated as the density difference between resident brine at the storage
AN US
layer top and storage layer bottom. The latter is determined by using, e.g., the formation thickness, 𝐻, and the geothermal gradient, 𝛻𝑇𝑔𝑒𝑜 , to derive the bottom temperature and pressure.
M
2.5.3 Viscosity submodel
ED
The dynamic viscosity model by Mao and Duan [59] was incorporated into the analysis tool to calculate the viscosity of NaCl solutions as a function of temperature, pressure and salinity.
PT
The viscosity model was chosen because of its applicability to the pressure and temperature
CE
conditions that are found under geological CO2 storage and a wide salinity range of 0-6 m NaCl. Viscosity is highly temperature dependent, decreasing with increasing temperature.
AC
Salinity also affects brine viscosity but to a lesser extent than temperature. The salinity effect on viscosity depends on the electrolyte; for NaCl solutions, the viscosity increases with ionic strength [59]. The influence of pressure is comparably small, as higher pressure raises the brine viscosity only slightly [60]. This study only addresses density-driven instability that is induced by ∆ρS and ∆ρT , and potential influences on instability that arise from viscosity changes that are caused by dissolved CO2 are not considered. Viscosity changes that follow CO2 dissolution are assumed to be negligible.
ACCEPTED MANUSCRIPT 2.5.4 Molecular diffusivity submodel The CO2 diffusivity, 𝐷, flattens concentration gradients within the porous medium. In this study, a dilute solution was assumed and the Wilke-Chang correlation [61], where the molecular diffusivity changes with temperature and solvent viscosity, 𝜇, was applied, as in
CR IP T
Eq. 4: ̃𝑤 )1/2 𝑇 7.4 ∙ 10−8 (𝜃𝑀 𝐷= 𝜇𝑉̅𝑆0.6
(4)
̃𝑤 is Here, 𝜃 is an empirical constant of the solvent, in this case water, with a value of 2.6. 𝑀
AN US
the molecular weight of water, and T is the temperature in Kelvin. The molar volume of the solute, i.e., CO2, 𝑉̅𝑆 , is 34.0 cm3/mol [61,62]. 2.5.5 Thermal diffusivity submodel
M
The thermal diffusivity, 𝜅, smooths temperature differences. 𝜅 is primarily temperature
ED
dependent, decreasing with temperature. The 𝜅 that is used to calculate 𝑅𝑎 is defined as the ratio between the thermal conductivity of the medium, 𝜆, and the fluid heat capacity, (𝑐𝜌)𝑓 ,
𝜅 = 𝜆 / (𝑐𝜌)𝑓
(5)
CE
PT
as in Eq. 5 [52]:
AC
Here, 𝑐𝑓 and 𝜌𝑓 are the heat capacity and density of the fluid, respectively. In this study, the heat capacity and density of water were calculated by employing the thermodynamic formulation of IAPWS97 [57], taking into account both pressure and temperature. The thermal conductivity, 𝜆, was estimated with the correction [63,64] in Eq. 6: 𝜆(𝑇) = 𝜆20°𝐶 − 10−3 ∙ (𝑇 − 293) ∙ (𝜆20°𝐶 − 1.38) ∙ [𝜆20°𝐶 ∙ (1.8 ∙ 10−3 𝑇)−0.25∙𝜆20°𝐶 + 1.28] ∙ 𝜆20°𝐶 −0.64
(6)
ACCEPTED MANUSCRIPT Here, the temperature, T, is given in Kelvin and the thermal conductivity at 20 °C, 𝜆20°𝐶 , is given in W/(m K). A sedimentary basin of sandstone with 𝜆20°𝐶 = 1.5 W/(m K) was assumed, a value in the range that was reported for sedimentary rocks [65]. Eq. 6 was found to be valid for both dry and liquid-saturated sandstones with porosities that are not too high [63]. In this study, liquid-saturated sandstone was considered.
CR IP T
2.5.6 Calculation of 𝑅𝑎, 𝑡𝐶 and 𝐹
The values of 𝜇, ∆𝜌𝑆 , ∆𝜌𝑇 , 𝐷 and 𝜅, which were calculated by the above submodels, were transferred (Fig. 2) to submodels to compute 𝑅𝑎, 𝑅𝑎𝑆 and 𝑅𝑎 𝑇 (Eq. 1), 𝑡𝐶 (Eq. 2) and 𝐹
AN US
(Table 2 and Eq. 3). All the scaling laws that are given in Table 2 are incorporated into the analysis tool.
2.6 Comparison study
M
The performance of each submodel that was described in sections 2.5.1 through 2.5.5 was
ED
verified against its original source study with excellent agreement. Furthermore, the analysis tool’s performance was compared to the study by Lindeberg and Wessel-Berg [8] as a
PT
benchmark by comparing calculated S-numbers for four different North Sea reservoir
CE
conditions. The S-number definition is given in the figure caption of Fig. 3. Because Lindeberg and Wessel-Berg [8] assumed a constant ∆𝜌𝑆 , a first comparison was made by
AC
replacing the temperature-, pressure- and salinity-dependent ∆𝜌𝑆 that was calculated by the analysis tool with their constant value. Fig. 3 shows that the calculated S-numbers concurred with the S-numbers that were estimated by Lindeberg and Wessel-Berg [8]. Only slight differences originated from the marginal differences in the thermal and molecular diffusivities and viscosities and the differences in the estimated ∆ρT . In the present study, lower values of ∆ρT were obtained as the analysis tool also accounts for the effect of the pressure gradient, which counteracts thermal expansion.
ACCEPTED MANUSCRIPT
To evaluate the effect of a non-constant ∆𝜌𝑆 -assumption, the pressure-, temperature- and salinity-dependent ∆𝜌𝑆 was calculated by the analysis tool and compared to the constant ∆𝜌𝑆 . The S-numbers that were obtained by the analysis tool showed similar trends as those that were obtained by Lindeberg and Wessel-Berg [8] but were systematically lower. This
CR IP T
indicates that the constant ∆𝜌𝑆 -assumption by Lindeberg and Wessel-Berg [8] overestimates S-numbers (case 1 by 21.5 %, case 2 by 40.1 %, case 3 by 52.2 % and case 4 by 57.4 %) for all conditions. The largest overestimations coincide with the highest temperature and pressure
AN US
conditions.
2.7 Sensitivity analysis
For this study, local and global sensitivity analyses were applied to identify influential
M
parameters on system variables. The analyses were conducted by connecting the developed
ED
analysis tool to iTOUGH2 [66] through the PEST protocol [68,69]. In the local sensitivity analysis, partial derivatives, i.e., the ratios between variable changes
PT
and the parameter perturbations that give rise to them, were scaled to obtain dimensionless
CE
sensitivity coefficients. The scaling was conducted by multiplying the partial derivatives with the ratio of the expected parameter variation and variable standard deviation. The
AC
dimensionless sensitivity coefficients permitted a comparison of the relative influence of parameters of different types, units or scales [66]. In addition, a global sensitivity analysis was conducted to investigate the parameter space, calculate sensitivity measures and identify parameter-interaction effects. The Morris one-at-a-time (OAT) method [67] was used, which partitions the parameter ranges into equally sized intervals and randomly generates sets of parameters. The generated parameter values are then perturbed and the changes in the variables are calculated. The analysis is performed over multiple “paths” (sets of parameter
ACCEPTED MANUSCRIPT values and orders of perturbing parameters) to obtain a global sensitivity index, the mean elementary effect (EE), for each parameter over the parameter space. A high absolute EE indicates a highly influential parameter, while a high standard deviation (SD) of the EE indicates the presence of non-linearity or interaction effects [66]. Furthermore, the uncertainty
CR IP T
in an EE estimate is given by 𝑆𝐷/𝑟 0.5, where 𝑟 is the number of paths [67].
2.8 Investigated cases
Analyses of the prerequisites for density-driven instability were conducted over broad ranges of conditions that are likely found under CO2 storage in deep saline aquifers, such as
AN US
temperatures of 30-130 °C, pressures of 7.85-60 MPa and aqueous NaCl solutions of 0-6 m (mol NaCl/kg H2O).
Trends in the CO2 solubility and subsequent ∆𝜌𝑆 over these broad ranges of conditions were
M
studied by using the analysis tool. Then, two aspects were studied to obtain an overview of
ED
the relative solutal and thermal contributions to density-driven instability over the ranges of these conditions. First, the trends of the relative solutal and thermal driving forces of
PT
convection, i.e., the buoyancy ratio, 𝑁 = ∆𝜌𝑆 /∆𝜌𝑇 , were studied within aquifers with different geothermal gradients, salinities and aquifer properties (Table 3). Second, the relative
CE
solutal and thermal contributions to the onset of instability, i.e., 𝑅𝑎, were studied through trends in 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 for the ranges of conditions that are given in Table 3. The ratio 𝜑𝑅𝑎𝑆 /
AC
𝑅𝑎 𝑇 is influenced by both 𝑁 and the diffusivity ratio. Thereafter, systematic calculations of 𝑅𝑎 over the broad thermodynamic and salinity ranges and permeabilities, which are given in Table 3, were conducted to obtain an overview of the 𝑅𝑎 range that is expected during CO2 storage in deep saline aquifers and to investigate the 𝑅𝑎 trends and limiting conditions for density-driven instability.
ACCEPTED MANUSCRIPT Furthermore, the impact of basin type on the prerequisites for density-driven instability was investigated; the 𝑅𝑎 and corresponding 𝑡𝐶 were calculated and analysed for the thermodynamic and salinity conditions that are likely found in a cold basin and a warm basin, respectively. Cold basins have low surface temperatures and geothermal gradients compared to warm basins, thus generally exhibiting lower temperatures at the same depths [70]. In
CR IP T
earlier studies [70,71], the basin type was shown to affect the CO2 storage capacity of free phase CO2; however, its impact on density-driven convective mixing, the focus of this study, has not been explored. The values of the surface temperature and geothermal gradients that are used in the present study (Table 3) were obtained from Bachu [71] and the 𝑡𝐶 by using Eq.
AN US
2.
To relate this study's theoretical findings to real cases, the 𝑅𝑎𝑆 and its corresponding 𝑡𝐶 for five CO2 pilot test and storage sites were calculated by using the parameters that are shown in
M
Table 4 and Eq. 2.
ED
Additionally, a procedure that permits the application of this study's results in a straightforward manner was developed to estimate 𝑅𝑎𝑆 for sites of one's own choosing while
PT
accounting for salinity and thermodynamic conditions. To create a basis for accurate and
CE
efficient site-specific estimations of 𝑅𝑎𝑆 , the thermodynamic- and salinity-dependent part of Eq. 1, i.e., ∆𝜌𝑆 /𝜇𝐷, hereafter referred to as the dynamic factor, was generated for a wide
AC
range of temperature (30-130 ˚C), pressure (10-60 MPa) and salinity (0-4 m NaCl) conditions. The impacts of the thermodynamic and salinity conditions and the basin type on the trends in 𝐹 and the duration of the convective regime were also studied. A comparison of 𝐹 that used both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws under cold basin conditions, which are given in Table 3, was conducted. The value of k that was used made the comparison of the scaling laws possible within their valid 𝑅𝑎𝑆 ranges, shown in Table 2. For the comparison, the
ACCEPTED MANUSCRIPT salinity was set to 0 m NaCl, which led to upper limit estimates of 𝐹. The impacts of the basin type and storage depth on the trends in 𝐹 were studied for three different salinities by using the scaling law by Slim [38], Eq. 3 and the properties in Table 3. The impacts of the basin type, salinity and storage depth on the duration of the convective regime were also
CR IP T
investigated by using the parameterization by Slim [38] in Table 2. Additionally, local and global sensitivity analyses were conducted to determine the relative influence of different field conditions and properties on density-driven instability and
convective mixing. In the local sensitivity analysis, the influences of the field parameters
AN US
(storage depth, basin type, i.e., ∇𝑇𝑔𝑒𝑜 , 𝜑, 𝑘, 𝐻 and salinity) on the variables 𝑥𝑆 , ∆𝜌𝑆 , 𝐷, 𝜇, 𝑅𝑎𝑆 , 𝐹, and 𝑡𝐶 and the duration of the convective regime were investigated. The analysis included 𝐹 that was based on both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws. However, only the influences of the field parameters on 𝑅𝑎𝑆 were analysed in the global
ED
analyses are given in Table 5.
M
sensitivity analysis (Morris OAT method). The parameters that were used in the sensitivity
PT
3. Results
CE
3.1 Trends in CO2 solubility and ∆𝜌𝑆 under broad ranges of conditions The CO2 solubility into brine and the subsequent ∆𝜌𝑆 % at different temperatures, pressures
AC
and salinities are shown in Fig. 4 and Fig. 5, respectively. Similar trends in CO2 solubility and ∆𝜌𝑆 % are produced (although the results deviate slightly at high temperatures and pressures at low salinities) by both solubility models [47,49]: decreasing trends with decreasing pressure and increasing salinity. Furthermore, the CO2 solubility first decreases with temperature until a pressure-dependent minimum is reached and then increases with temperature, while ∆𝜌𝑆 % decreases with increasing temperature. However, the impact of the temperature and pressure
ACCEPTED MANUSCRIPT conditions on both CO2 solubility and ∆𝜌𝑆 % diminishes with increasing salinity. Because of the salinity and temperature validity limits of the solubility models, only results from the Spycher and Pruess [49] or Duan and Sun [47] solubility model are shown for salinities over 4 m NaCl or temperatures over 100 °C, respectively. Fig. 5 shows that ∆𝜌𝑆 % is generally around 1 % to a few tenths of a percentage for the wide range of conditions that were studied.
CR IP T
The maximum ∆𝜌𝑆 %, ~2 % (not shown), corresponds to thermodynamic conditions of 30 °C, 60 MPa and 0 m NaCl. However, at high temperatures and salinities, ∆𝜌𝑆 % is negative, i.e., a density decrease follows CO2 dissolution into brine. This result is consistent with the result by Lu et al. [74] on CO2 storage efficiency, who showed the temperature, or “equal density
AN US
temperature”, beyond which CO2-saturated brine is lighter than resident brine (∆𝜌𝑆 < 0) to be 118 °C at 20 MPa and 4 m salinity. However, for the wider salinity range up to 6 m NaCl that was investigated in the present study, Fig. 5 shows that the “equal density temperature”
M
decreases with increasing salinity for a certain pressure; at conditions of 20 MPa and 6 m
ED
NaCl, the “equal-density temperature” is 80 °C and is even lower at lower pressures.
instability
PT
3.2 The relative solutal and direct thermal contributions to the onset of
CE
Previous studies [8,46,51] have found that the thermal contribution to convection is negligible compared to the solutal contribution; ∆𝜌𝑇 is generally one order of magnitude less than ∆𝜌𝑆
AC
[51] and 𝑅𝑎 𝑇 is orders of magnitude less than 𝑅𝑎𝑆 [8]. In the present study, a relative comparison of ∆𝜌𝑆 and ∆𝜌𝑇 , i.e., 𝑁, was made within an aquifer for different ∇𝑇𝑔𝑒𝑜 and salinities (Fig. 6). The assumed 𝐻 in this study (100 m) is within the range of 𝐻 that is considered in previous studies [8,51]. The 𝑁 varies widely, depending on the thermodynamic and salinity conditions and increases as the ∇𝑇𝑔𝑒𝑜 and salinity decrease. The highest ratios are obtained for thermodynamic conditions that correspond to shallow depths. Over the wide
ACCEPTED MANUSCRIPT ranges of conditions that were studied, ∆𝜌𝑇 varies from one order of magnitude smaller than ∆𝜌𝑆 to reaching equal magnitude or even larger. At larger depths, the ratios, which vary with ∇𝑇𝑔𝑒𝑜 , turn negative as a result of the above-mentioned induced density decrease. Fig. 7 shows 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 for the same ∇𝑇𝑔𝑒𝑜 and salinities that were used to study 𝑁. The trend in φRaS /RaT with depth is a combined effect of the trends in both N and the diffusivity ratio.
CR IP T
These trends are similar, but the influence of the latter is greater. 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 increases with decreasing ∇𝑇𝑔𝑒𝑜 , depth, 𝜑 and salinity. The ratio is generally large, with 𝑅𝑎𝑆 being several hundreds to tens of thousand times larger than 𝑅𝑎 𝑇 . However, the ratio turns negative at a
AN US
certain depth, which is determined by ∇𝑇𝑔𝑒𝑜 , and high salinities.
3.3 Trends in 𝑅𝑎 under broad ranges of conditions
In Fig. 8, 𝑅𝑎 values for different temperature, pressure, salinity and permeability conditions
M
are shown. The Y-axis scale of each graph is different. The displayed trends in 𝑅𝑎 are of greater importance than the absolute values because the absolute values are based on likely,
ED
but assumed, formation properties such as 𝐻 and 𝜑. The trends in 𝑅𝑎 are similar to those that
PT
are displayed by ∆𝜌𝑆 % in Fig. 5: a decreasing trend with increasing temperature and salinity and an increasing trend with increasing pressure. The pressure influence on 𝑅𝑎 diminishes
CE
with temperature at high salinities (≥4 m NaCl). Negative values of 𝑅𝑎 correspond to negative
AC
∆𝜌𝑆 %, which is mentioned in connection to Fig. 5. The solid red lines in Fig. 8 show the 𝑅𝑎𝐶 range 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 . Values of 𝑅𝑎 that fall above the upper limit of the 𝑅𝑎𝐶 range indicate prerequisites for the onset of instability under the specific thermodynamic and salinity conditions. The 𝑘 has a large influence on the magnitude of 𝑅𝑎. In relatively high 𝑘 (10-12 m2) and low-moderate salinity aquifers, 𝑅𝑎 is generally large, on the order of a thousand to tens of thousands, i.e., 𝑅𝑎 ≫ 𝑅𝑎𝐶 , and the prerequisites for the onset of densitydriven instability are good. For even higher 𝑘 or thicker layers, 𝑅𝑎 values of a few hundreds
ACCEPTED MANUSCRIPT of thousands can be expected. However, even in aquifers with high 𝑘, 𝑅𝑎 values that are close to zero or negative can be found at high salinity and certain thermodynamic conditions. For aquifers with lower k (10-14 m2) and low-moderate salinity, 𝑅𝑎 is one order higher or of the same order of magnitude as the RaC range. Under these conditions, instability is highly dependent on the prevailing thermodynamic and salinity conditions. In low 𝑘 (10-15 m2)
CR IP T
aquifers, the 𝑅𝑎 values are generally small, falling below or within the 𝑅𝑎𝐶 range, and the prerequisites for the onset of density-driven instability are limited. Furthermore, Fig. 8 shows that 𝑅𝑎 values in high salinity aquifers are generally low, except at high 𝑘 and low-moderate temperatures. In addition, large negative values of 𝑅𝑎 are obtained for high salinity and
AN US
temperature conditions at high 𝑘.
3.5 Trends in 𝑅𝑎 and 𝑡𝐶 under field conditions - impact of the basin type (i.e.,
M
∇𝑇𝑔𝑒𝑜 )
ED
In Fig. 9, 𝑅𝑎 values for different 𝑘 and salinities for a cold basin at 800-3500 m depth and a warm basin at 800-2400 m depth are shown. For the warm basin case, results were obtained
PT
down to a depth of 1800 m by using the Spycher solubility model [49] and at lower depths down to 2400 m by using the Duan and Sun solubility model [47] because of the temperature
CE
validity limits of the solubility models that were used. Furthermore, negative 𝑅𝑎 values are
AC
not shown because of the logarithmic Y-axis scales. A general trend is that the highest 𝑅𝑎 is obtained for temperature and pressure conditions that correspond to shallow depths of around 800 m, and then 𝑅𝑎 decreases with depth. Furthermore, the 𝑅𝑎 values are generally higher in the cold basin than in the warm basin because of the smaller ∇𝑇𝑔𝑒𝑜 of the former, which induces lower temperatures and thus larger ∆𝜌𝑆 (Fig. 5) and leads to larger 𝑅𝑎 in turn. In the cold basin, 𝑅𝑎 generally exceeds the
ACCEPTED MANUSCRIPT assumed 𝑅𝑎𝐶 range at high 𝑘 (10-13 m2 to 10-12 m2), indicating prerequisites for the onset of instability, except at high salinities. At high salinities, the depth determines the prerequisites for density-driven instability. Similar trends are seen in the warm basin, but the depths at which the prerequisites for instability exist are shallower (at even lower salinities) than in the cold basin. For aquifers in the intermediate 𝑘 range (10-14 m2 to 10-13 m2) with low salinities
CR IP T
(≤1 m NaCl) and for both basin types, the 𝑅𝑎 values generally exceed the 𝑅𝑎𝐶 range for the studied depths. However, as the salinity increases, the influence of depth on the onset of
instability increases. In lower 𝑘 (10-15 m2 to 10-14 m2) aquifers, the prerequisites for densitydriven instability are more limited and are highly affected by the salinity, depth and basin type
AN US
(i.e., ∇𝑇𝑔𝑒𝑜 ). For formations with low 𝑘 (≤ 10-15 m2), the 𝑅𝑎 values are small. Only for low salinity and shallow depths in the cold basin do the 𝑅𝑎 values fall within the considered 𝑅𝑎𝐶 range but do not exceed it, i.e., the prerequisites for instability are very limited for low 𝑘.
M
In Fig. 10, the impact of salinity on 𝑅𝑎𝑆 is shown under cold and warm basin conditions. The
ED
𝑅𝑎𝑆 values that are obtained under cold basin conditions generally exceed the 𝑅𝑎𝐶 range for low-moderate salinities. However, for high salinities, the depth, i.e., thermodynamic
PT
conditions, highly influences whether prerequisites for instability exist. For a cold basin with a salinity of ~5 m NaCl, the temperature and pressure conditions for which 𝑅𝑎𝑆 first falls
CE
within the 𝑅𝑎𝐶 range correspond to a depth of ~3300 m but decreases to 2600 m for a salinity
AC
of ~6 m NaCl. Under warm basin conditions, the depths at which 𝑅𝑎𝑆 first falls within the 𝑅𝑎𝐶 range are even shallower than those under cold basin conditions and occur at even lower salinities. In the warm basin, only salinities <3 m NaCl give rise to values of 𝑅𝑎𝑆 that exceed the 𝑅𝑎𝐶 range within the whole depth range. As in the cold basin, the depth has a major influence on the prerequisites for the onset of instability at higher salinities. For the warm basin and salinities of ~4 m NaCl, ~5 m NaCl and ~6 m NaCl, 𝑅𝑎𝑆 first falls within the 𝑅𝑎𝐶 range at a depth of approximately 1800 m, 1400 m and 1100 m, respectively.
ACCEPTED MANUSCRIPT In Fig. 11, the 𝑡𝐶 values for different salinities under cold and warm basin conditions are shown. These values of 𝑡𝐶 are based on the values of 𝑅𝑎𝑆 in Fig. 10. Only 𝑡𝐶 for 𝑅𝑎𝑆 ≥13.8 are shown. Furthermore, the values of 𝑡𝐶 were calculated by using Eq. 2 and the numerical constant 𝑐0 =31, which represents impermeable upper boundary conditions [30], thus giving the minimum estimates of 𝑡𝐶 compared to assuming a permeable upper boundary condition.
CR IP T
One trend (Fig. 11) is that 𝑡𝐶 increases with salinity and depth and the influence of the depth increases with salinity. However, 𝑡𝐶 is generally smaller for cold basin conditions than for warm basin conditions at the same depth.
AN US
3.6 Field sites and a new procedure for the accurate site-specific estimation of 𝑅𝑎𝑆
The calculated values of 𝑅𝑎𝑆 and 𝑡𝐶 for five CO2 pilot test and storage sites, both off- and
M
onshore, are given in Table 6. These pilot test and storage sites have different temperature, pressure and salinity settings and formation properties, as given in Table 4. The largest 𝑅𝑎𝑆 is
ED
obtained for the Sleipner site; the ∆𝜌𝑆 % is almost 1 %, considerably higher than for the other
PT
sites. The 𝑅𝑎𝑆 for Frio, Ketzin and Snøhvit are of the same order of magnitude. The lowest 𝑅𝑎𝑆 is obtained for the In Salah site in the Krechba field. This 𝑅𝑎𝑆 , which falls below the 𝑅𝑎𝐶
CE
range, occurred because of the combined effect of small values of 𝑘, 𝐻, and ∆𝜌𝑆 under the thermodynamic and salinity conditions of the site. The 𝑡𝐶 values for the sites are given in
AC
Table 6. The difference in the 𝑡𝐶 between the Ketzin and Snøhvit sites, despite the sites’ similar 𝑅𝑎𝑆 and 𝐻 values, arises from the difference in 𝐷. A two-step procedure was developed to apply this study's results in a straight-forward manner to estimate 𝑅𝑎𝑆 for a site of one's own choosing: 1) Given the temperature, pressure and salinity of a specific site, Fig. 12 can be used to obtain the dynamic factor, i.e., ∆𝜌𝑆 /𝜇𝐷 in Eq. 1, for a wide range of conditions (proxy-equations for the dynamic factor can be found in
ACCEPTED MANUSCRIPT Appendix A), and approximations of the dynamic factor for salinities that are not shown can be made by interpolation between two dynamic factors that are obtained for a lower and higher salinity. 2) Once the dynamic factor is obtained, the site-specific 𝑅𝑎𝑆 (Eq. 1) can be estimated by dividing by 𝜑 and multiplying with 𝐻 (in unit m), 𝑘 (in unit m2) and 𝑔 (in unit m/s2), which varies with layer inclination. The estimated 𝑅𝑎𝑆 can then serve as an indicator
CR IP T
for the prerequisite for the onset of density-driven instability at a site and used to evaluate the suitability of storage sites.
3.7 Trends in 𝐹 under field conditions and impacts on the duration of the
AN US
convective regime
In Fig. 13, the 𝐹 values that are predicted by different scaling law types (𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent) and assumed upper boundary conditions (Table 2) are shown under cold
M
basin conditions. Both the non-dimensional mass flux, 𝐹 ∗ , and corresponding dimensional mass flux, 𝐹, are shown. The mass fluxes are obtained for zero salinity conditions and
ED
therefore constitute an upper bound estimate of the flux for the given thermodynamic
PT
conditions. Fig. 13 gives an overview of the implications of using a specific type of scaling law or upper boundary condition on the predicted 𝐹. The comparison shows consistency in
CE
the predicted trend in 𝐹 whether 𝑅𝑎𝑆 -dependent or 𝑅𝑎𝑆 -independent scaling laws were used; 𝐹 increased with increasing depth. In Fig. 13, the experimentally derived 𝑅𝑎𝑆 -dependent
AC
scaling laws are shown with their uncertainty intervals (dotted lines). For the range of 𝑅𝑎𝑆 values that were studied, the steady convective flux that was predicted by the 𝑅𝑎𝑆 independent scaling laws for impermeable (by, e.g., Slim [38] and Elenius et al. [30]) and partially (0.2) permeable upper boundary fall within the uncertainty interval of the 𝑅𝑎𝑆 dependent scaling law that was experimentally derived by Neufeld et al. [32]. A relatively wide range of predicted flux sizes is shown in Fig. 13 for different upper permeability
ACCEPTED MANUSCRIPT boundary conditions. The smallest steady convective fluxes result from the 𝑅𝑎𝑆 -dependent scaling law that was experimentally derived by Backhaus et al. [33]. Results are only shown for cold basin conditions because similar results were obtained for warm basin conditions. The effects of salinity, basin type and storage depth on the trends in 𝐹 (using the scaling law by Slim [38]) down to a depth of 2500 m and the corresponding ∆𝜌𝑆 /𝜇 ratio are shown in Fig.
CR IP T
14. Only values of 𝐹 at conditions for which the system’s 𝑅𝑎𝑆 is within the valid 𝑅𝑎𝑆 range (Table 2) are shown in Fig. 14, which is the reason why 𝐹 under warm basin conditions (at 3 m NaCl) is not shown for the entire depth range. A substantial decrease in 𝐹 is visible as
AN US
salinity increases in both cold and warm basins, but the basin type (i.e., ∇𝑇𝑔𝑒𝑜 ) also has a visible impact on 𝐹. Furthermore, the trends in 𝐹 closely follow the salinity-dependent trends in the ∆𝜌𝑆 /𝜇 ratio. At low salinities (0-1 m NaCl), 𝐹 and ∆𝜌𝑆 /𝜇 increase with depth. However, as salinity increases, the trend in 𝐹 with depth also changes. The trend in 𝐹 evolves,
M
first increasing with depth, then becoming more uniform, and finally decreasing with depth as
ED
salinity increases.
In addition, Fig. 15 shows the effects of salinity, basin type and storage depth on the
PT
convective flux regime with regards to 𝐹 and the 𝑅𝑎𝑆 -dependent duration (when applying the
CE
parameterization by Slim [38]). Only the convective flux regime is shown, not the preceding or successive flux regimes. Fig. 15 shows that combinations of factors such as different
AC
salinity, storage depth and basin type (i.e., ∇𝑇𝑔𝑒𝑜 ) highly affect 𝐹 and the approximate onset time and duration of the convective regime; e.g., high salinity, large storage depth and warm basin conditions (large ∇𝑇𝑔𝑒𝑜 ), which correspond to a low 𝑅𝑎𝑆 , delay the onset of convection, decrease the magnitude of the steady convective flux and shorten the duration of the convective regime.
3.8 The relative influence of different field conditions and properties
ACCEPTED MANUSCRIPT The local sensitivity analysis, where scaled sensitivity coefficients were used to compare the influence of different field conditions and properties on the 𝑥𝑆 , ∆𝜌𝑆 , 𝐷, 𝜇, 𝑅𝑎𝑆 , 𝐹, 𝑡𝐶 and duration of the convective regime, has provided several findings (Fig. 16). The analysis indicates that 1) 𝑅𝑎𝑆 is sensitive to all the studied parameters. However, 2) 𝑘, 𝜑 and 𝐻 have the highest influence on 𝑅𝑎𝑆 . 3) The 𝑥𝑆 , ∆𝜌𝑆 , 𝐷 and 𝜇 are sensitive to ∇𝑇𝑔𝑒𝑜 , salinity and
CR IP T
storage depth. 4) The most influential parameters on 𝐹 when applying both 𝑅𝑎𝑆 -independent (𝐹 a ) and 𝑅𝑎𝑆 -dependent (𝐹 b ) scaling laws are 𝑘 and salinity. However, 5) for 𝐹 a, the
sensitivity to 𝜑 and 𝐻 differs, which is expected because 𝑅𝑎𝑆 is highly sensitive to these parameters. 6) The ∇𝑇𝑔𝑒𝑜 and storage depth seem to have an influence on 𝐹. 7) The most
AN US
influential parameters on 𝑡𝐶 are 𝑘 and 𝜑, followed by salinity, ∇𝑇𝑔𝑒𝑜 and storage depth. 8) All the studied parameters influence the duration of the convective regime, but the most influential are 𝑘, 𝐻, 𝜑 and salinity.
M
The result from the global sensitivity analysis (Morris OAT method) regarding the influence
ED
of the field conditions and properties on 𝑅𝑎𝑆 is shown in Fig. 17. The analysis shows that 𝑘, salinity, 𝜑 and 𝐻 have the highest influence on 𝑅𝑎𝑆 , as indicated by the high absolute values
PT
of their mean EE. However, all the studied parameters display a statistically significant non-
CE
zero impact on 𝑅𝑎𝑆 as they fall below the dashed lines that mark the standard error of the mean, 𝑆𝐷/𝑟 0.5 (see figure caption). Furthermore, high values of SD EE indicate the presence
AC
of non-linearity or interaction effects, and the high influence of salinity on 𝑅𝑎𝑆 , although not visible in the local sensitivity analysis, becomes evident when the average parameter influence over the entire parameter range is considered.
4. Discussion When investigating the prerequisites for enhanced solubility trapping during geological CO2 storage under broad ranges of conditions, we found that the driving force for convective
ACCEPTED MANUSCRIPT mixing ceases under certain thermodynamic and salinity conditions; e.g., under certain conditions, CO2 dissolution into brine causes a density decrease instead of an increase, which concurs with the findings of Lu et al. [74]. Even if the density of the CO2-saturated brine is lighter than the underlying resident brine, this value is still heavier than that of the overlying CO2-rich phase, which leads to a stable density gradient in the formation as long as the CO2-
point, a risk for upward buoyancy-driven flow may exist.
CR IP T
rich phase remains in the formation. However, if lighter CO2-saturated brine reaches a spill
In investigating the relative solutal and direct thermal contributions to the onset of instability
AN US
under broad ranges of conditions, ∆𝜌𝑇 and 𝑅𝑎 𝑇 were found to be highly dependent on the ∇𝑇𝑔𝑒𝑜 and 𝐻 but were generally smaller than ∆𝜌𝑆 and 𝑅𝑎𝑆 , except under certain thermodynamic conditions at high salinity (Fig. 6 and Fig. 7). However, the absolute magnitudes of 𝑅𝑎𝑆 and 𝑅𝑎 𝑇 at these high salinity conditions are both small, resulting in a low
M
𝑅𝑎 that falls slightly above, within or below the 𝑅𝑎𝐶 range. The direct thermal contribution
ED
during these conditions can be quantified by the change in critical depth (depth at which 𝑅𝑎 = 4𝜋 2 ) between when 𝑅𝑎 𝑇 is considered in Eq. 1 and when it is not. Accounting for the direct
PT
thermal contribution leads to an increase in critical depth. For the cold and warm basin conditions (6 m NaCl, 𝜑=0.3, 𝐻=100 m) that are shown in Fig. 9, the critical depths increase
CE
by ~44 m and ~54 m, respectively. Thus, the direct thermal impact, i.e., the increase in the
AC
storage depth for which a prerequisite for density-driven instability exists, is limited to at most a few tens of meters. Several previous studies [8,46,51] have concluded that the thermal contribution to convection is negligible compared to that of the solutal. However, we found that clear, non-negligible, indirect thermal contributions to the onset of instability and convective mixing exist according to the results of the conducted trend and sensitivity analyses, although the direct thermal contribution generally is small. The trends in 𝑅𝑎 (Fig. 9), 𝑅𝑎𝑆 (Fig. 10), 𝑡𝐶 (Fig. 11) and 𝐹 (Fig. 14) were shown to be affected by the basin type
ACCEPTED MANUSCRIPT (i.e., the ∇𝑇𝑔𝑒𝑜 ). The local sensitivity analysis (Fig. 16) identified that the ∇𝑇𝑔𝑒𝑜 influences the 𝑥𝑆 , 𝜇, 𝐷, ∆𝜌𝑆 , 𝑅𝑎𝑆 , 𝐹, 𝑡𝐶 and the convective regime duration. In addition, the global sensitivity analysis confirmed that the ∇𝑇𝑔𝑒𝑜 has a statistically significant non-zero impact on 𝑅𝑎𝑆 .
CR IP T
Investigating the prerequisites for density-driven instability under broad ranges of field conditions (Fig. 8) indicated that combinations of high temperature (≥ 80 °C), moderate-high salinity (≥ 4 m NaCl) and low 𝑘 (≤ 10-14 m2) disfavour density-driven instability. 𝑘 has a
major influence on the magnitude of 𝑅𝑎, which increases with 𝑘. This was also confirmed by
AN US
the local and global sensitivity analyses (Figs. 16 and 17) in this study. However, this study also shows that the trends in 𝑅𝑎 over the broad ranges of thermodynamic and salinity conditions follow the trends in ∆𝜌𝑆 , namely, decreasing with temperature and salinity and increasing with pressure. Thus, deep saline aquifers with relatively low temperature and
M
salinity and high 𝑘 provide the optimal conditions for the onset of density-driven instability,
ED
such as the conditions that were reported for the North Sea reservoir [8]. Indeed, the result from the site-specific 𝑅𝑎𝑆 comparison showed that the Sleipner site has exceptionally good
PT
prerequisites for density-driven instability compared to the other sites that were studied.
CE
Furthermore, cold basin conditions are more favourable to the onset of density-driven instability than warm basin conditions, rendering higher 𝑅𝑎 (Fig. 9) and shorter 𝑡𝐶 (Fig. 11)
AC
for the same depth. Thus, the lowest depth at which the prerequisites for instability exist is shallower for warm basin conditions than for cold basin conditions (Figs. 9 and 10), given the same salinity. Consequently, as this depth decreases with increasing salinity, the choice of storage depth is especially important for CO2 storage in warm basins with high salinity. The general trend that was found in this study for all basin types, salinities, and 𝑘 is that the largest 𝑅𝑎 is obtained at temperature and pressure conditions that correspond to relatively
ACCEPTED MANUSCRIPT shallow depths around 800 m and thereafter decreases with depth. At the same time, the 𝑡𝐶 was found to increase with storage depth for moderate-high salinities (≥3 m NaCl), but this trend ceases with decreasing salinity (Fig. 11), which indicates that CO2 storage at shallower depths favours the onset of density-driven instability and short 𝑡𝐶 , especially within highly saline aquifers. In the study by Bachu [71], the optimal storage depths for free phase CO2
CR IP T
based on capacity and drilling costs were concluded to be 800-1000 m for cold basins and 1500-2000 m for warm basins. For cold basins, the suggested depth range also coincides with the largest 𝑅𝑎. Accordingly, this depth range is optimal from both a capacity (free phase
storage) and onset of instability (enhanced solubility trapping) perspective. For warm basins,
AN US
the optimal depth that was mentioned by Bachu [71] does not coincide with the depth for the largest 𝑅𝑎. The values of 𝑅𝑎 that are found in this depth range only exceed the 𝑅𝑎𝐶 range for salinities <4 m NaCl. These results indicate that convection might not set in at this depth
M
range at higher salinities, instead favouring the onset of instability and short 𝑡𝐶 to store CO2 at shallow depths. However, these shallow depth conditions do not necessarily create the highest
ED
steady convective flux, 𝐹. As shown by the result from the local sensitivity analysis (Fig. 16)
PT
and trend analyses, 𝐹 is highly dependent on the 𝑘 (comparing Figs. 13 and 14), the salinity (Figs. 14 and 16), the amount of exchanges with the capillary transition zone (the upper
CE
boundary permeability, Fig. 13), the basin type (i.e., ∇𝑇𝑔𝑒𝑜 ) and the storage depth (Fig. 14). Both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws show that 𝐹 follows the salinity-
AC
dependent trend of the ∆𝜌𝑆 /𝜇 ratio with depth within an aquifer. However, at moderate to high salinities (≥3 m NaCl), the thermodynamic conditions that benefit both 𝐹 and 𝑡𝐶 coincide at shallow depth if the prerequisites for instability exist. In addition, the results (Figs. 10 and 14) imply that the storage sites with the highest 𝑅𝑎𝑆 do not necessarily have the highest 𝐹.
ACCEPTED MANUSCRIPT The impact of salinity, basin type and storage depth on the size of 𝐹 and the approximate onset and duration of the steady convective regime is compelling (Fig. 15), i.e., the time scale of solubility trapping is substantially affected by the prevailing thermodynamic and salinity conditions of the storage site and is thus of great concern to risk assessment during site characterization.
CR IP T
The findings of this study can be applied in a straight-forward manner to estimate the 𝑅𝑎𝑆 of a specific site or to obtain a greater understanding of which field conditions promote densitydriven instability and convection and enhance solubility trapping. Site-specific estimates of 𝑅𝑎𝑆 can be obtained by using the procedure that was presented in section 3.6 (applying the
AN US
proxy-equations in Appendix A or the result in Fig. 12) without running the developed analysis tool, making it possible to easily evaluate and compare storage locations’
prerequisites for the onset of instability, i.e., enhanced solubility trapping, as part of a site
M
characterization effort. The findings of this study could also be used for future studies together with global field data to map storage formations with conditions that are favourable
ED
for enhanced solubility trapping.
PT
The present study considers density-driven instability and convection in a homogenous isotropic porous layer, not convection that is driven by changes in, e.g., 𝜇 from CO2
CE
dissolution or temperature differences. The direct thermal contribution to convection can be
AC
expected to increase when thermal effects on 𝜇 are also considered. The effect of dissolved CO2 on 𝜇 in aqueous NaCl solutions was previously studied [75,76] and a viscosity model was presented [77]. Neglecting the aqueous CO2 concentration dependency of 𝜇 will result in a slight underestimation of the aqueous phase 𝜇 (increasing with aqueous CO2 concentration), a disregard of an indirect effect on 𝐷 and a slight overestimation of the 𝑅𝑎𝑆 . This was also confirmed by calculations of RaS (results not shown) that were conducted with the abovementioned viscosity model incorporated into the analysis tool. However, including the effect
ACCEPTED MANUSCRIPT of aqueous CO2 concentration on 𝜇 only slightly affects the values of 𝑅𝑎𝑆 , while the trend in 𝑅𝑎𝑆 remains the same, which confirms that the influence of thermodynamic and salinity conditions on 𝜇 are stronger than that of the aqueous CO2 concentration. Furthermore, NaCl was selected to represent the salinity in this study because Na+ is one of
CR IP T
the most commonly cations that are found in groundwater [78]. However, other ions may be present in the formation fluid or dominate at certain sites. These ions influence CO2 solubility; both the ionic strength and the composition (valence and size) of the ions determine the
salting-out effect. Thus, brine that contains divalent cations (e.g., MgCl2 or CaCl2 solutions)
AN US
instead of monovalent cations (e.g., NaCl or KCl solutions) exhibits a lower CO2 solubility at the same concentration [79]. The properties of the aqueous phase, e.g., 𝜇, are also affected by ions. The mentioned effects are expected to influence density-driven convection. Further research could aim to quantify these effects over broad ranges of field conditions and the
M
effect of co-injected impurities on density-driven convection, which can be expected to affect,
PT
5. Conclusions
ED
e.g., the solubility and thereby the induced density changes driving convection.
This study provides an overview of the field conditions for which the prerequisites for
CE
density-driven instability and convection, which enhance solubility trapping, exist. An analysis was conducted over broad thermodynamic, salinity and permeability conditions that
AC
are relevant for CO2 storage in deep saline aquifers. New insight into the trends of 𝑅𝑎, 𝑡𝐶 and 𝐹 over these broad ranges of field conditions is provided and the key parameters that control these are identified. For this purpose, a new analysis tool for comprehensive calculations of 𝑅𝑎, 𝑡𝐶 and 𝐹 was developed to account for the temperature, pressure and salinity dependency of ∆𝜌𝑆 , ∆𝜌𝑇 , 𝜇, 𝐷 and 𝜅 and their combined effect on density-driven instability. These dependences are usually disregarded or simplified but, as this study shows, are of great
ACCEPTED MANUSCRIPT importance to understand the trends of 𝑅𝑎, 𝑡𝐶 and 𝐹 under field conditions. In addition, the relative influence of field parameters, e.g. storage depth and basin type (i.e., geothermal gradient), were analysed through local and global sensitivity analyses. In our study of the conditions that favour and disfavour density-driven instability and convective mixing, we found that the magnitude of 𝑅𝑎 is highly influenced by permeability
CR IP T
but that the trends of 𝑅𝑎 follow the trends of ∆𝜌𝑆 , i.e., decreasing with temperature and salinity and increasing with pressure.
The direct thermal contribution (∆𝜌𝑇 ) to density-driven instability was also found to be small
AN US
compared to the solutal contribution (∆𝜌𝑆 ), except at high salinities. Contrary to previous studies, this study shows that the geothermal gradient has a non-negligible effect on densitydriven instability and convective mixing when accounting for both direct and indirect thermal effects. Cold basin conditions generally render higher 𝑅𝑎 compared to warm. However, the
M
largest 𝑅𝑎 for all basin types and salinities is obtained for thermodynamic conditions that
ED
correspond to relatively shallow depths around 800 m, which then decreases with depth. This indicates that shallow depth CO2 storage, especially at high salinities, favours the onset of
PT
density-driven instability and short 𝑡𝐶 . Shallow depth conditions do not, however, necessarily
CE
render the largest 𝐹, as shown in a trend analysis of both 𝑅𝑎𝑆 -dependent and 𝑅𝑎𝑆 -independent scaling laws for 𝐹; instead, the optimal storage depth for large 𝐹 depends on the salinity,
AC
owing to the ratio ∆𝜌𝑆 /𝜇. Altogether, the salinity, basin type and storage depth highly affect the onset of density-driven instability, 𝑡𝐶 , 𝐹 and the duration of the convective regime, subsequently altering the solubility trapping time scale. In addition, a straight-forward and efficient procedure to estimate site-specific 𝑅𝑎𝑆 was presented. 𝑅𝑎𝑆 may be obtained for any site by means of the dynamic factor (∆𝜌𝑆 /𝜇𝐷), which
ACCEPTED MANUSCRIPT is comprehensively calculated by the developed analysis tool. The estimated 𝑅𝑎𝑆 can serve as an indicator for density-driven instability prerequisites at a site when evaluating its suitability as a storage location.
CR IP T
Acknowledgments The authors would like to thank Shide Mao at the China University of Geosciences for
providing the FORTRAN source code to his viscosity model [59] and Nicholas Spycher at Lawrence Berkeley National Laboratory for providing the FORTRAN source code to his
AN US
solubility model [48,49]. The research that led to these results has received funding from the European Community's 7th Framework Programme MUSTANG and PANACEA projects (grant agreement numbers 227286 and 282900) and Swedish Research Council (VR). The
M
authors would like to thank the anonymous reviewers for their reviews and constructive
ED
comments, which improved this manuscript. Appendix A
PT
Proxy-equations for the dynamic factors are shown in Fig. 12.
CE
Dynamic factor = (𝑎1 + 𝑎2 ∙ T 2 + 𝑎3 ∙ T 3 + 𝑎4 ∙ T 4 + 𝑎5 ∙ 𝑃 + 𝑎6 ∙ P2 + 𝑎7 ∙ P3 + 𝑎8 ∙ P4 + 𝑎9 ∙ T ∙ P + 𝑎10 ∙ 𝑇 2 ∙ 𝑃 + 𝑎11 ∙ 𝑇 ∙ 𝑃2 ) ∙ 1012
Here, the pressure, P, is in bars and the temperature, T, is in ˚C.
AC
The coefficients are given in Table A.
ACCEPTED MANUSCRIPT References [1] IPCC. IPCC Special Report on carbon dioxide capture and storage. Cambridge University Press, Cambridge (2005). [2] Garcia JE. Density of aqueous solutions of CO2. LBNL-49023. Lawrence Berkeley
CR IP T
National Laboratory, Berkeley (2001). [3] Pruess K, Zhang K. Numerical modeling studies of the dissolution-diffusion-convection process during CO2 storage in saline aquifers. LBNL-1243E. Lawrence Berkeley National
AN US
Laboratory, Berkeley (2008).
[4] Ennis-King J, Paterson L. Rate of dissolution due to convective mixing in the underground storage of carbon dioxide. 6th International Conference on Greenhouse Gas Control
M
Technologies 1: 507-510 (2003).
[5] Yang C, Gu Y. Accelerated mass transfer of CO2 in reservoir brine due to density-driven
PT
2430-2436 (2006).
ED
natural convection at high pressures and elevated temperatures. Ind Eng Chem Res 45(8):
[6] Hassanzadeh H, Pooladi-Darvish M, Keith DW. Modelling of convective mixing in CO2
CE
storage. J Can Petrol Technol 44(10): 43-50 (2005).
AC
[7] Weir GJ, White SP, Kissling WM. Reservoir storage and containment of greenhouse gases. Transp Porous Med 23: 37-60 (1996). [8] Lindeberg E, Wessel-Berg D. Vertical convection in an aquifer column under a gas cap of CO2. Energ Convers Manage 38: Suppl., S229-S234 (1997).
ACCEPTED MANUSCRIPT [9] Farajzadeh R, Salimi H, Zitha PLJ, Bruining H. Numerical simulation of density-driven natural convection in porous media with application for CO2 injection projects. Int J Heat Mass Tran 50: 5054-5064 (2007). [10] Kneafsey TJ, Pruess K. Laboratory flow experiments for visualizing carbon dioxide-
CR IP T
induced, density-driven brine convection. Transp Porous Med 82: 123-139 (2010). [11] Tilton N, Riaz A. Nonlinear stability of gravitationally unstable, transient, diffusive boundary layers in porous media. J Fluid Mech 745: 251-278 (2014).
[12] Elenius MT, Johannsen K. On the time scales of nonlinear instability in miscible
AN US
displacement porous media flow. Comput Geosci 16: 901-911 (2012).
[13] Bestehorn M, Firoozabadi, A. Effect of fluctuations on the onset of density-driven
M
convection in porous media. Phys Fluids 24: 114102 (2012).
[14] Riaz A, Hesse M, Tchelepi HA, Orr Jr. FM. Onset of convection in a gravitationally
ED
unstable diffusive boundary layer in porous media. J Fluid Mech 548: 87-111 (2006).
PT
[15] Hassanzadeh H, Pooladi-Darvish M, Keith DW. Stability of a fluid in a horizontal saturated porous layer: effect of non-linear concentration profile, initial, and boundary
CE
conditions. Transp Porous Med 65: 193-211 (2006).
AC
[16] Rapaka S, Chen S, Pawar RJ, Stauffer PH, Zhang D. Non-modal growth of perturbations in density-driven convection in porous media. J Fluid Mech 609: 285-303 (2008). [17] Xu X, Chen S, Zhang D. Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv Water Resour 29: 397-407 (2006). [18] Ennis-King J, Paterson L. Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations. SPE J, September 2005: 349- 356 (2005).
ACCEPTED MANUSCRIPT [19] Ennis-King J, Preston I, Paterson L. Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys Fluids 17: 084107 (2005). [20] Rapaka S, Pawar RJ, Stauffer PH, Zhang D, Chen S. Onset of convection over a transient base-state in anisotropic and layered porous media. J Fluid Mech 641: 227-244 (2009).
driven convection. Transp Porous Med 72: 241-253 (2008).
CR IP T
[21] Hong JS, Kim MC. Effect of anisotropy of porous media on the onset of buoyancy-
[22] Hassanzadeh H, Pooladi-Darvish M, Keith DW. Scaling behavior of convective mixing,
AN US
with application to geological storage of CO2. AIChE J 53(5): 1121-1131 (2007).
[23] Azin R, Raad SMJ, Osfouri S, Fatehi R. Onset of instability in CO2 sequestration into saline aquifer: scaling relationship and the effect of perturbed boundary. Heat Mass Transfer
M
49: 1603-1612 (2013).
[24] Green CP, Ennis-King J. Effect of vertical heterogeneity on long-term migration of CO2
ED
in saline formations. Transp Porous Med 82: 31-47 (2010).
PT
[25] Pau GSH, Bell JB, Pruess K, Almgren AS, Lijewski MJ, Zhang K. High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers. Adv
CE
Water Resour 33(4): 443-455 (2010).
AC
[26] Chen C, Zeng L, Shi L. Continuum-scale convective mixing in geological CO2 sequestration in anisotropic and heterogeneous saline aquifers. Adv Water Resour 53: 175-187 (2013).
[27] Ennis-King J, Paterson L. Coupling of geochemical reactions and convective mixing in the long-term geological storage of carbon dioxide. Int J Greenhouse Gas Control 1: 86-93 (2007).
ACCEPTED MANUSCRIPT [28] Javaheri M, Abedi J, Hassanzadeh H. Onset of convection in CO2 sequestration in deep inclined saline aquifers. J Can Petrol Technol 48(8): 22-27 (2009). [29] Slim AC, Ramakrishnan TS. Onset and cessation of time-dependent, dissolution-driven convection in porous media. Phys Fluids 22: 124103. DOI:10.1063/1.3528009 (2010).
CR IP T
[30] Elenius MT, Nordbotten, JM, Kalisch, H. Effects of a capillary transition zone on the stability of a diffusive boundary layer. IMA J Appl Math 77: 771-787 (2012).
[31] Emami-Meybodi H, Hassanzadeh H. Two-phase convective mixing under a buoyant
AN US
plume of CO2 in deep saline aquifers. Adv Water Resour 76: 55-71 (2015).
[32] Neufeld JA, Hesse MA, Riaz A, Hallworth MA, Tchelepi HA, Huppert HE. Convective dissolution of carbon dioxide in saline aquifers. Geophys Res Lett 37: L22404 (2010).
M
[33] Backhaus S, Turitsyn K, Ecke RE. Convective instability and mass transport of diffusion layers in a Hele-Shaw geometry. Phys Rev Lett 106(10): 104501 (2011).
ED
[34] Zhang W, Li Y, Omambia AN. Reactive transport modeling of effects of convective
PT
mixing on long-term CO2 geological storage in deep saline formations. Int J Greenhouse Gas Control 5: 241-256 (2011).
CE
[35] Hidalgo JJ, MacMinn CW, Juanes R. Dynamics of convective dissolution from a
AC
migrating current of carbon dioxide. Adv Water Resour 62: 511-519 (2013). [36] Szulczewski ML, Hesse MA, Juanes R. Carbon dioxide dissolution in structural and stratigraphic traps. J Fluid Mech 736: 287-315 (2013). [37] Slim AC, Bandi MM, Miller JC, Mahadevan L. Dissolution-driven convection in a HeleShaw cell. Phys Fluids 25: 024101 (2013).
ACCEPTED MANUSCRIPT [38] Slim AC. Solutal-convection regimes in a two-dimesional porous medium. J Fluid Mech 741: 461-491 (2014). [39] Hesse MA. Mathematical modeling and multiscale simulation for CO2 storage in saline aquifers. PhD, Dept. of Energy Resources Engineering, Stanford University (2008).
media. Phys Rev Lett 109: 264503 (2012).
CR IP T
[40] Hidalgo JJ, Fe J, Cueto-Felgueroso L, Juanes R. Scaling of convective mixing in porous
[41] Ghesmat K, Hassanzadeh H, Abedi J. The impact of geochemistry on convective mixing in a gravitationally unstable diffusive boundary layer in porous media: CO2 storage in saline
AN US
aquifers. J Fluid Mech 673: 480-512 (2011).
[42] Islam A, Korrani AKN, Sepehrnoori K, Patzek T. Effects of geochemical reaction on
Mass Tran 77: 519-528 (2014).
M
double diffusive natural convection of CO2 in brine saturated geothermal reservoir. Int J Heat
ED
[43] Farajzadeh R, Ranganathan P, Zitha PLJ, Bruining J. The effect of heterogeneity on the character of density-driven natural convection of CO2 overlying a brine layer. Adv Water
PT
Resour 34: 327-339 (2011).
CE
[44] Ranganathan P, Farajzadeh R, Bruining H, Zitha PLJ. Numerical simulation of natural convection in heterogeneous porous media for CO2 geological storage. Transp Porous Med
AC
95(1): 25-54 (2012).
[45] Kong X-Z, Saar MO. Numerical study of the effects of permeability heterogeneity on density-driven convective mixing during CO2 dissolution storage. Int J Greenhouse Gas Control 19: 160-173 (2013).
ACCEPTED MANUSCRIPT [46] Islam AW, Lashgari HR, Sephernoori K. Double diffusive natural convection of CO2 in a brine saturated geothermal reservoir: study of non-modal growth of perturbations and heterogeneity effects. Geothermics 51: 325-336 (2014). [47] Duan Z, Sun R. An improved model calculating CO2 solubility in pure water and aqueous
(2003).
CR IP T
NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chem Geol 193(3-4): 257-271
[48] Spycher N, Pruess K, Ennis-King J. CO2-H2O mixtures in the geological sequestration of CO2. I. Assessment and calculation of mutual solibilities from 12 to 100˚C and up to 600 bar.
AN US
Geochim Cosmochim Acta 67(16): 3015-3031 (2003).
[49] Spycher N, Pruess K. CO2-H2O mixtures in the geological sequestration of CO2. II. Partitioning in chloride brines at 12-100˚C and up to 600 bar. Geochim Cosmochim Acta
M
69(13): 3309-3320 (2005).
ED
[50] Loodts V, Rongy L, De Wit A. Impact of pressure, salt concentration, and temperature on the convective dissolution of carbon dioxide in aqueous solutions. CHAOS 24: 043120
PT
(2014).
CE
[51] Javaheri M, Abedi L, Hassanzadeh H. Linear stability analysis of double-diffusive convection in porous media, with application to geological storage of CO2. Transp Porous
AC
Med 84: 441-456 (2010). [52] Nield DA, Bejan A. Convection in porous media, 3rd edition. Springer (2006). [53] Kumar A, Ozah R, Noh M, Pope GA, Bryant S, Sepehrnoori K, et al., Reservoir simulation of CO2 storage in deep saline aquifers. SPE J 10(3): 336-348 (2005).
ACCEPTED MANUSCRIPT [54] Li D, Graupner BJ, Bauer S. A method for calculating the liquid density for the CO2H2O-NaCl system under CO2 storage condition. Energy Procedia 4: 3817-3824 (2011). [55] Mao S, Duan Z. The P,V,T,x properties of binary aqueous chloride solutions up to T=573 K and 100 MPa. J Chem Thermodyn 40(7): 1046-1063 (2008).
100 MPa. Energ Fuel 22(3): 1666-1674 (2008).
CR IP T
[56] Duan Z, Hu J, Li D, Mao S. Densities of the CO2-H2O-NaCl systems up to 647 K and
[57] Wagner W, Cooper JR, Dittmann A, Kijima J, Kretzschmar HJ, Kruse A, et al., The IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. J
AN US
Eng Gas Turbines Power 122(1): 150-182 (2000).
[58] Archer DG, Wang P. The dielectric constant of water and Debye-Hückel limiting law
M
slopes. J Phys Chem Ref Data 19(2): 371-411 (1990).
[59] Mao S, Duan Z. The viscosity of aqueous alkali-chloride solutions up to 623 K, 1,000
ED
bar, and high ionic strength. Int J Thermophys 30(5): 1510-1523 (2009).
PT
[60] Adams JJ, Bachu S. Equations of state for basin geofluids: algorithm review and intercomparison for brines. Geofluids 2(4): 257-271 (2002).
CE
[61] Wilke CR, Chang P. Correlation of diffusion coefficients in dilute solutions. AIChE J
AC
1(2): 264-270 (1955).
[62] Cussler EL. Diffusion: mass transfer in fluid systems. 3rd edition, Cambridge University Press, Cambridge (2009). [63] Somerton WH. Thermal properties and temperature-related behavior of rock/fluid systems. Developments in petroleum science, 37. Elsevier, Amsterdam (1992).
ACCEPTED MANUSCRIPT [64] Tikhomirov V. Conductivity of rocks and their relationship with density, saturation and temperature. Neft Khoz 46(4): 36 (1968). [65] Lee Y, Deming D. Evaluation of thermal conductivity temperature corrections applied in terrestrial heat flow studies. J Geophys Res 103(B2): 2447-2454 (1998).
CR IP T
[66] Finsterle S. iTOUGH2 User’s Guide. LBNL-40040. Lawrence Berkeley National Laboratory, Berkeley (2007).
[67] Morris M D. Factorial sampling plans for preliminary computational experiments.
AN US
Technometrics 33(2): 161-174 (1991)
[68] Doherty J, PEST: Model-independent parameter estimation. Watermark numerical computing, Brisbane, Australia (2008), http://pesthomepage.org/
M
[69] Finsterle S, Zhang Y. Solving iTOUGH2 simulation and optimization problems using the PEST protocol. Environmental Modelling & Software 26: 959-968 (2011).
ED
[70] Bachu S. Sequestration of CO2 in geological media in response to climate change: road
PT
map for site selection using the transform of the geological space into the CO2 phase space. Energ Convers Manage 43(1): 87-102 (2002).
CE
[71] Bachu S. Screening and ranking of sedimentary basins for sequestration of CO2 in
AC
geological media in response to climate change. Environ Geol 44(3): 277-289 (2003). [72] Michael K, Allinson G, Golab A, Sharma S, Shulakova V. CO2 storage in saline aquifers II - experience from existing storage operations. Energy Procedia. 1: 1973-1980 (2009). [73] Trémosa J, Castillo C, Quang Vong C, Kervévan C, Lassin A, Audigane P. Long-term assessment of geochemical reactivity of CO2 storage in highly saline aquifers: Application to Ketzin, In Salah and Snøhvit storage sites. Int J Greenhouse Gas Control 20: 2-26 (2014).
ACCEPTED MANUSCRIPT [74] Lu C, Han WS, Lee S-Y, McPherson BJ, Lichtner PC. Effects of density and mutual solubility of a CO2-brine system on CO2 storage in geological formations: "Warm" vs. "cold" formations. Adv Water Resour 32(12): 1685-1702 (2009). [75] Bando S, Takemura F, Nishio M, Hihara E, Akai M. Viscosity of aqueous NaCl solutions
(2004).
CR IP T
with dissolved CO2 at (30 to 60) °C and (10 to 20) MPa. J Chem Eng Data 49: 1328-1332
[76] Fleury M, Deschamps H. Electrical conductivity and viscosity of aqueous NaCl solutions with dissolved CO2. J Chem Eng Data 53: 2505-2509 (2008).
AN US
[77] Islam AW, Carlson ES. Viscosity models and effects of dissolved CO2. Energ Fuel 26: 5330-5336 (2012).
[78] Freeze RA, Cherry JA. Groundwater. Upper Saddle River, Prentice Hall, Inc. (1979). [79] Zerai B, Saylor B Z, Allen DE. Quantification of CO2 trapping and storage capacity in
M
the subsurface: uncertainty due to solubility models. In: McPherson BJ, Sundquist ET,
ED
editors. Carbon sequestration and its role in the global carbon cycle. Geophysical Monograph
AC
CE
PT
Series 183, American Geophysical Union. p. 249-260 (2009).
ACCEPTED MANUSCRIPT
M
AN US
CR IP T
Figure captions
Fig. 1. Schematic of a homogenous, isotropic porous layer. The symbols and abbreviations are
AC
CE
PT
ED
defined in the nomenclature.
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 2. Schematic of the developed analysis tool and its overall structure and submodels.
AN US
CR IP T
ACCEPTED MANUSCRIPT
∆𝜌𝑆 𝐷
+
∆𝜌𝑇 𝑔𝑘𝐻
M
Fig. 3. Comparison of S-numbers, 𝑆 = (
𝜅
)
𝜇
, from the study by Lindeberg and
Wessel-Berg [8] and those that were calculated by the developed analysis tool, which account
AC
CE
PT
ED
for constant and varying ∆𝜌𝑆 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 4. Solubility of CO2, 𝑥𝑆 , at different temperatures and pressures for salinities of (a) 1 m
M
NaCl, (b) 4 m NaCl and (c) 6 m NaCl. The solid and dashed lines illustrate the results of the
AC
CE
PT
ED
solubility models by Duan and Sun [47] and Spycher and Pruess [49], respectively.
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 5. Density changes from CO2 dissolution, ∆𝜌𝑆 %, at different temperatures and pressures
M
for salinities of (a) 1 m NaCl, (b) 4 m NaCl and (c) 6 m NaCl. The solid and dashed lines
AC
CE
PT
[49], respectively.
ED
illustrate the results of the solubility models by Duan and Sun [47] and Spycher and Pruess
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 6. Buoyancy ratios, 𝑁 = ∆𝜌𝑆 /∆𝜌𝑇 , at various depths for different geothermal gradients at
AC
CE
PT
ED
M
salinities of (a) 0 m NaCl and (b) 6 m NaCl. Hydrostatic pressure is assumed.
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 7. The ratio 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 at various depths for different geothermal gradients at salinities
AC
CE
PT
ED
M
of (a) 0 m NaCl and (b) 6 m NaCl. Hydrostatic pressure is assumed.
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
Fig. 8. 𝑅𝑎 for different temperature, pressure, salinity and permeability conditions. The left column shows the results for k=10-12 m2, the middle column for k=10-14 m2 and the right column for k=10-15 m2. The coloured curves indicate different pressures: black - 7.85 MPa, green - 30 MPa, and blue - 60 MPa. The solid and dashed lines illustrate the results of the solubility models by Duan and Sun [48] and Spycher and Pruess [49], respectively. The solid red lines indicate the 𝑅𝑎𝐶 range: 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 9. 𝑅𝑎 for cold basin (black lines) and warm basin (magenta lines) conditions with ∇𝑇𝑔𝑒𝑜 = 25 °C/km and 50 °C/km, respectively. 𝑇𝑠𝑢𝑟𝑓 =10 °C and hydrostatic pressure is assumed. 𝑅𝑎
ED
values are shown for different permeabilities and salinities of 0 m NaCl, 1 m NaCl, 4 m NaCl
AC
CE
PT
and 6 m NaCl. The solid red lines indicate the 𝑅𝑎𝐶 range: 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 10. 𝑅𝑎𝑆 at different salinities and various depths. To the left: for a cold basin with ∇𝑇𝑔𝑒𝑜 = 25 °C/km. To the right: for a warm basin with ∇𝑇𝑔𝑒𝑜 = 50 °C/km. For both cases, 𝑇𝑠𝑢𝑟𝑓 =10 °C
ED
and k=10-13 m2, and hydrostatic pressure is assumed. The solid red lines indicate the 𝑅𝑎𝐶
AC
CE
PT
range: 13.8 < 𝑅𝑎𝐶 < 4𝜋 2 .
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 11. The 𝑡𝐶 at different salinities and various depths. To the left: for a cold basin with
AC
CE
PT
ED
𝑇𝑠𝑢𝑟𝑓 =10 °C and k=10-13 m2.
M
∇𝑇𝑔𝑒𝑜 =25 °C/km. To the right: for a warm basin with ∇𝑇𝑔𝑒𝑜 =50 °C/km. For both cases,
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
CE
Fig. 12. Dynamic factor, ∆𝜌𝑆 /𝜇𝐷𝐶𝑂2 , for different temperatures, pressures and salinities. The
AC
arrow indicates the direction of increasing pressure. Pressure increases with 10 MPa for each curve. The first curve corresponds to 10 MPa, the second curve corresponds to 20 MPa and so on up to 60 MPa.
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 13. Comparison of the scaling laws of the steady convective flux. To the left is the non-
M
dimensional mass transfer (𝐹 ∗ ) at zero salinity for a cold basin with ∇𝑇𝑔𝑒𝑜 =25 °C/km,
ED
𝑇𝑠𝑢𝑟𝑓 =10 °C, k=10-12 m2 and a hydrostatic pressure gradient. The dotted lines represent
AC
CE
PT
experimental uncertainty bounds. To the right is the corresponding dimensional flux (𝐹).
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 14. To the left is the dimensional steady convective flux, 𝐹, at three different salinities and various depths within a cold basin (black symbols) and a warm basin (magenta symbols)
ED
with ∇𝑇𝑔𝑒𝑜 =25 °C/km and 50 °C/km, respectively. 𝑇𝑠𝑢𝑟𝑓 =10 °C and k=10-13 m2, and
AC
CE
PT
hydrostatic pressure is assumed. To the right is the corresponding ∆𝜌𝑆 /𝜇.
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 15. Convective regime duration at two different depths and salinities for a cold basin (solid lines) and a warm basin (dashed lines) with ∇𝑇𝑔𝑒𝑜 =25 °C/km and 50 °C/km,
AC
CE
PT
ED
respectively, and a hydrostatic pressure gradient. 𝑇𝑠𝑢𝑟𝑓 =10 °C and k=10-13 m2.
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 16. Scaled sensitivity coefficients, which show the relative influence of different field parameters on 𝑥𝑆 , ∆𝜌𝑆 , 𝐷, 𝜇, 𝑅𝑎𝑆 , 𝐹, 𝑡𝐶 and the convective regime duration. 𝐹 a is based on the
AC
CE
PT
Neufeld et al. [32].
ED
𝑅𝑎𝑆 -independent scaling law by Slim [38]. 𝐹 b is based on the 𝑅𝑎𝑆 -dependent scaling law by
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 17. Mean EE, which shows the relative influence of different field parameters on 𝑅𝑎𝑆 .
M
SD EE conveys non-linearity or interaction effects. The dashed lines represent the mean EE
AC
CE
PT
ED
±2 (SD/𝑟 0.5).
ACCEPTED MANUSCRIPT
Table captions Table 1. 𝑅𝑎𝐶 for a homogenous and isotropic porous medium under different upper boundary conditions. For all cases, the lower boundary is assumed to be impermeable to fluid and heat/solute flow, with a prescribed temperature and/or solute concentration. 𝑅𝑎𝐶 4𝜋 27.10
impermeable to fluid and heat/solute flow, constant 𝑇 and 𝐶𝑆 .a
26.09
partially permeable (0.1) to fluid flowb, constant 𝐶𝑆 .c
23.45
partially permeable (0.5) to fluid flowb, constant 𝐶𝑆 .c
21.76
fully permeable to fluid flowb, constant 𝐶𝑆 .c
impermeable to fluid flow, constant 𝑃, 𝑇 and 𝐶𝑆 .a
CR IP T
Upper boundary conditions
2
AN US
13.8 impermeable to fluid flow, constant 𝑃 and 𝐶𝑆 .d a Based on linear stability analysis for thermally and solutally induced instability with a steady base state (initial temperature and concentration distributions that vary linearly with depth) [52]. b
Continuity of flux and pressure. The relative permeability of the capillary transition zone is given within brackets. c
Numerically derived for solutally (CO2) induced instability with an evolving base state [29].
d
M
Minimum value that is obtained from (energy) stability boundary mapping for solutally (CO2) induced instability with an evolving base state [29,37].
F*
Valid 𝑅𝑎𝑆 range
Num.
impermeable
-
Num.
[25]
impermeable
-
Num.[38]
impermeable
102-5∙104
Num.[30]
impermeable
Method
0.017
CE
0.017 to 0.018 0.017 a 0.02
AC
0.025 0.044
Upper boundary
[39]
PT
convective flux, 𝐹 ∗ .
ED
Table 2. Numerically and experimentally derived scaling laws of the non-dimensional steady
0.075 (0.12 ± 0.03)𝑅𝑎𝑆 −0.16±0.02
Num.
[38]
10 -5∙104
partially permeable (0.6)c
102-5∙104
partially permeable (0.2)
Num.[38] Num.
[30]
Analogue fluid exp.
c
fully permeable [32]
-
2
b
3
10 -106
(0.045 ± 0.025)𝑅𝑎𝑆 −0.24±0.06 Analogue fluid exp.[33] 6∙103-9∙104 a ∗ Duration of the convective regime is approximated by the parameterization 𝐹 = 0.017, 1100 < 𝑡 < 16 ∙ 𝑅𝑎𝑆 by Slim [38], where 𝑡 is the non-dimensional time that is obtained by multiplying the dimensional time with (𝑔𝑘∆𝜌𝑆 )2 (𝜇𝜑)2 𝐷
.
b
Full exchange with the capillary transition zone.
c
Relative permeability of the capillary transition zone.
ACCEPTED MANUSCRIPT
Table 3. Aquifer properties and conditions that were used to calculate 𝑁, 𝜑𝑅𝑎𝑆 /𝑅𝑎 𝑇 , 𝑅𝑎 and 𝐹 for cold and warm basins. 𝑇𝑠𝑢𝑟𝑓 =10 °C, 𝜑=0.3 and 𝐻=100 m.
𝑁 and 𝜑Ra S /RaT
1.0∙10
a
𝛻𝑇𝑔𝑒𝑜 [°C/km]
Salinity [m NaCl]
-
25, 30, 40, 50
0-6 0-6
-15
1.0∙10 -1.0∙10
-12
7.85-60
30-130
50
-15
-12
a
-
25, 50
1.0∙10 -1.0∙10
𝐹-scaling laws 𝐹 trends
T [°C]
a
1.0∙10-12
a
-13
a
1.0∙10 A hydrostatic pressure gradient is assumed.
0-6
25
0
25, 50
0, 1, 3
AN US
Ra Basin types
P [MPa]
-12
CR IP T
𝑘 [m2]
Cases
Table 4. Properties of the CO2 pilot test and storage sites that were used as input parameters to calculate the site-specific 𝑅𝑎𝑆 and 𝑡𝐶 .
𝑘 [m ] 2
Ketzin73
24
80
1.5∙10
𝜑 [-]
0.3
T [°C]
56
P [MPa]
-12
15.2
a
Sleipner73
M
H [m]
Frio73
7.5∙10
-13
0.23
In Salah74
Snøhvit74
20
75
250 5.0∙10
-12
1.0∙10
-14
8.83∙10-13
0.37
0.15
0.15
34
37
96
98
7.3
10.3
18.0
28.5
ED
Parameters
CE
PT
Salinity [m NaCl] 1.598 4.288 0.587 2.712 2.566 a Salinities converted to m NaCl from TDS, taking into account the site-specific thermodynamic conditions.
AC
Table 5. Parameters that were used in the sensitivity analyses. 𝑘 [m2]
Depth [m]
𝜑 [-]
1.0∙10-13
1000
0.3
𝐻 [m]
∇𝑇𝑔𝑒𝑜 [°C/km]
Salinity [m NaCl]
25
1
a
Local analysis initial value
100
Global analysisb range 1.0∙10-15-1.0∙10-11 800-1800 0.1-0.4 10-100 25-50 0-6 a The initial parameter values were perturbed by 10 % and centred finite differences were used. The sensitivity coefficients were scaled with the ratio of the parameter and variable variations, which were assumed to be 10 % of their initial values. b
Parameter ranges were partitioned into 11 equal intervals and 800 paths were used (= 5600 simulations).
ACCEPTED MANUSCRIPT Table 6. Calculated ∆𝜌𝑆 and 𝑅𝑎𝑆 for specific CO2 pilot test and storage sites.
∆𝜌𝑆 [kg/m3] ∆𝜌𝑆 % [%] 𝑅𝑎𝑆 [-] 𝑡𝐶 [days] a
Frio
Ketzin
Sleipner
In Salah
Snøhvit
5.3
2.7
10.1
1.5
2.0
0.51
0.24
0.99
0.14
0.19
3113.9
3748.0
177,702.5
8.6
3870.9
0.3
c
6.3
102.0
23.8
c 𝑡𝐶 b [days] 29.5 480.5 1.3 a Corresponding to permeable upper boundary condition, 𝑐0 = 31, [30].
CR IP T
112.2
b
Corresponding to impermeable upper boundary condition, 𝑐0 = 146, [14].
c
𝑡𝐶 is not calculated since stable conditions prevail, i.e., 𝑅𝑎𝑆 < 𝑅𝑎𝐶 .
AN US
Table A. Coefficients for the proxy-equations. 0.0 7.948136110915495 -0.003441435234038 0.000040935045132 -0.000000139741972 0.018651779590719 -0.000047204955860 0.000000078802468 -0.000000000046039 0.000006178398246 -0.000000060718778 -0.000000034826392 0.99820
0.5 6.707809604348991 -0.002930182594837 0.000034867612722 -0.000000118953377 0.015363104245509 -0.000037208693130 0.000000062337161 -0.000000000036472 0.000000757544439 -0.000000067646659 -0.000000026725659 0.99832
1.0 5.670286672377516 -0.002501054900298 0.000029768688332 -0.000000101467421 0.012661286761162 -0.000029034768917 0.000000048870428 -0.000000000028640 -0.000003567264417 -0.000000071268675 -0.000000020138373 0.99845
m NaCl 𝑎1 𝑎2 𝑎3 𝑎4 𝑎5 𝑎6 𝑎7 𝑎8 𝑎9 𝑎10 𝑎11 𝑅2
1.5 4.801002897984335 -0.002140719463234 0.000025483216074 -0.000000086761195 0.010429841104933 -0.000022323006341 0.000000037804293 -0.000000000022195 -0.000006990943279 -0.000000073206522 -0.000000014734211 0.99858
2.0 4.071765097530080 -0.001837906559726 0.000021879189254 -0.000000074386526 0.008580958668878 -0.000016808787401 0.000000028705856 -0.000000000016889 -0.000009651589099 -0.000000074281849 -0.000000010321768 0.99871
2.5 3.459177590320337 -0.001583147388553 0.000018845297093 -0.000000063964377 0.007044357126234 -0.000012279544781 0.000000021228394 -0.000000000012520 -0.000011663239757 -0.000000075013892 -0.000000006753049 0.99883
3.0
3.5
4.0
ED
PT
CE
AC m NaCl
M
m NaCl 𝑎1 𝑎2 𝑎3 𝑎4 𝑎5 𝑎6 𝑎7 𝑎8 𝑎9 𝑎10 𝑎11 𝑅2
ACCEPTED MANUSCRIPT
2.509486329629760 -0.001187411168486 0.000014129731257 -0.000000047755012 0.004691213083351 -0.000005513538486 0.000000010056460 -0.000000000005972 -0.000014106958036 -0.000000076689295 -0.000000001686914 0.99907
2.142827424385575 -0.001034325676134 0.000012304981938 -0.000000041478565 0.003790733063294 -0.000003017413526 0.000000005938987 -0.000000000003547 -0.000014685146127 -0.000000078002192 -0.000000000009779 0.99917
CR IP T
2.943800700123073 -0.001368518645321 0.000016288158757 -0.000000055176110 0.005763182036722 -0.000008561884080 0.000000015089251 -0.000000000008926 -0.000013121673473 -0.000000075742015 -0.000000003907715 0.99895
AC
CE
PT
ED
M
AN US
𝑎1 𝑎2 𝑎3 𝑎4 𝑎5 𝑎6 𝑎7 𝑎8 𝑎9 𝑎10 𝑎11 𝑅2