Amir Taheri ^{1,*}, Ole Torsæter ^{1}, Erik Lindeberg ^{2}, Nanji J. Hadia ^{3}, Dag WesselBerg ^{4}

Department of Geoscience and Petroleum, NTNU, Trondheim, Norway

SINTEF Petroleum Research, Trondheim, Norway

Institute of Chemical and Engineering Science (ICES), Singapore

Department of Mathematical Sciences, NTNU, Trondheim, Norway
Academic Editor: Sina R. Gomari
Special Issue: Applied Geological CO_{2} Storage and Geochemistry
Received: October 19, 2020  Accepted: February 07, 2021  Published: February 25, 2021
Advances in Chemical Research 2021, Volume 3, Issue 1, doi:10.21926/acr.2101012
Recommended citation: Taheri A, Torsæter O, Lindeberg E, Hadia NJ, WesselBerg D. Effect of Convective Mixing Process on Storage of CO_{2} in Saline Aquifers with Layered Permeability. Advances in Chemical Research 2021;3(1):21; doi:10.21926/acr.2101012.
© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.
Abstract
Convective mixing of freephase CO_{2} and brine in saline aquifers is an established technique to accelerate the CO_{2} dissolution process. Correct estimation of the convection onset time and rate of CO_{2} dissolution into brine are two crucial parameters regarding safety issues, as the timescale for dissolution corresponds to the same time over which the freephase CO_{2} has a chance to leak out from the storage site. In real practice, underground formations are heterogeneous with a layered structure, but the convective mixing in heterogeneous porous media has received less attention than the homogeneous one. This study aims to develop a basic understanding of the role of layered permeability media (layered structure with variation in permeability vertically) on the behavior of convective mixing via wellcontrolled laboratory experiments. The effects of layering and layer properties on the rate of dissolution of CO_{2} in water and geometries of the formed convection fingers are studied using a precise experimental setup with layeredpermeability HeleShaw cell geometry. Qualitative (snapshots of convection fingers) and quantitative data (amount of the dissolved CO_{2} into water) are collected simultaneously for a better understanding of the process. The behavior of convection fingers (after the onset of convection) and the effects of model properties on this mixing process are also discussed.
Keywords
Saline aquifer; layered permeability; convection; dissolution; HeleShaw cell
1. Introduction
Carbon dioxide (CO_{2}) storage in saline aquifers of huge volume has been recognized as an effective option for decreasing CO_{2} emission to the atmosphere. However, leakage of CO_{2} from these storage sites is one of the main concerns. The associated risk to this leakage is reduced by understanding the trapping mechanisms of CO_{2} into the brine and by taking suitable measures. One of these trapping mechanisms is the dissolution of supercritical CO_{2} into underground formation water, which is considered as a medium to longterm trapping mechanism. Over a long period of time, the injected CO_{2,} which has formed a thin layer of freephase CO_{2} below the caprock, gradually diffuses and dissolves into brine. The CO_{2} dissolution into brine, subsequently causes an increase in density of the brineCO_{2} solution, which overlies the less dense brine. This situation initiates gravitational instability and eventually leads to densitydriven natural convection and increasing the dissolution rate of free phase CO_{2} into the brine [1,2,3]. Densitydriven natural convection occurs when the Rayleigh dimensionless number $R_a=(\Delta\rho gkh)⁄(\varphi\mu D)$ is larger than 40, approximately [1]. The timing of the onset of this instability and the dissolution rate across the phase contact of CO_{2} and brine solution are the two crucial operational issues when assessing the feasibility of a potential storage site.
Numerical studies for describing the convective mixing process and its effectiveness in storing CO_{2} into the homogeneous saline aquifers are reported several years ago [2,3,4,5,6,7,8,9,10,11,12]. Furthermore, the gravitational instability of a diffusive boundary layer has been described using theoretical methodologies in recent years to predict the onset time of convection and the related unstable wavelength in the homogeneous models [13,14,15,16,17,18,19,20,21,22].
A limited number of experimental studies exist on the effectiveness of the convective mixing process in the homogeneous media regarding the storage of CO_{2} in geological formations [2,6,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. The occurrence of natural convection and the effectiveness of the convection mechanism (on the dissolution of CO_{2 }into brine) are investigated by recording the change in pressure in a bulk media, using a cylindrical PVT cell [25,26,27,28,33], and in a homogeneous porous media, using cylindrical core [6,30,31], where a fixed volume of CO_{2} is overlaying a column of distilled water. Convective currents and dynamics of convective mixing are visualized at ambient conditions using Schlieren and pH indicator techniques, for bulk modules of gasliquid systems, like CO_{2}water [23,24,33], and for homogeneous HeleShaw cells without porous media [2,30,31,32,34,35] and with porous media [29,30,31,32,35] by using analog or CO_{2}water fluids. The initiation, formation, and growth of small convective fingers and progress of this process as a function of time are discussed extensively [36]. The essential parameters of the convective mixing process, i.e., onset time of convection, critical wavelength, shape and growth of fingers over time, and the amount of dissolved CO_{2} into water, are measured qualitatively (by snapshots of convective mixing of CO_{2} and water) and quantitatively (by amounts of dissolved CO_{2} into water) in different homogeneous HeleShaw cell geometries with Rayleigh numbers varying in the range of 7009500 [36].
While heterogeneity and layered structure are the main characteristics of real geological formations, a limited number of studies exist about the gravitational instability of a diffusive boundary layer in such systems. Similar to homogeneous models, the onset conditions of natural convection in heterogeneous media (with simple geometry and boundary conditions) are also studied using linear stability analysis [37,38]. Effect of various types of heterogeneity, e.g., regular and irregular barrier patterns in porous media or multilayered porous media, on convective mixing are analyzed using numerical methods. In heterogeneous models with the distribution of horizontal impermeable barriers, the onset of convection is found to occur much sooner than the equivalent anisotropic homogeneous model. However, the constant dissolution flux (in the constant flux regime and after the onset of the instability) is not sensitive to barrier properties and depends only on the effective average properties [3,8,10,35,39]. It is also reported that the rate of CO_{2} dissolution in these types of heterogeneous media is higher than the equivalent homogeneous media, and the numerical simulations in equivalent homogenous porous media often underestimate the mass transfer rate of CO_{2} into water [11]. In multilayered porous media, though a toplayer with high permeability is more favorable for fast mass transfer, the natural convection process occurs with a toplayer of low permeability [40]. In such models with horizontal parallel layers and vertical variation of permeability, if the top layer has higher permeability, then the top layer properties govern the convective mixing in both strong and weak heterogeneity models. However, it is also seen that the toplayer with lower permeability can also be the governing one in the models with weak heterogeneity. Interestingly it is observed that in a model with a strong heterogeneity, the top layer with lower permeability has more effect on the amount of CO_{2} dissolved in brine than the bottom layer, and the impact of the bottom layer with lower permeability is increased with an increase in the thickness (bottom layer) [3]. Considering the behavior of convective mixing process, multilayered models are converged to homogeneous anisotropic models by increasing the number of layers [41]. The existence of low permeability rectangular structures inside a high permeability body results in different free convective processes at different spatial and temporal scales [42].
Experimental studies on convective mixing of CO_{2} and brine in twolayered porous media by quantifying pressure decline in a column of gas on top of a porous medium show that existence of a high permeable top layer is more desirable. This configuration causes an immediate start of convective mixing, but its strength decreases when the fingers penetrate into the low permeability layer [40]. It is also examined that strengthening the microheterogeneity in porous media accelerates the CO_{2} dissolution [42] and the existence of low permeability layers (embedded between relatively higher permeability zones) causes significant diffusive mixing in the low permeability layers than the convective mixing [44].
Despite numerous analytical and numerical studies on the convective mixing of CO_{2} and brine, insufficient experimental studies exist on the investigation of the accelerated mass transfer rate of CO_{2} into saline aquifers, especially in heterogeneous porous media with layered permeability. This work aims to develop a basic understanding of the mechanism of convective flow in layered permeability models, using experimental tools for identifying the behavior of convection fingers that meet the interface of the layers. For this purpose, the previous experimental setup of convective mixing studies [36,43] is upgraded using twolayered HeleShaw cells, and different tests are performed using CO_{2} and water at atmospheric conditions. Recording qualitative (snapshots of convective mixing) and quantitative measurements (pressure, temperature, and pump volumes) simultaneously give a better understanding. The timelapse movies are made from the snapshots that are collected during the experiments, which are suitable for improving public awareness. Lastly, the experimental results are compared qualitatively with numerical simulation models performed by Eclipse100.
2. Experimental Setup and Procedure
2.1 Experimental Setup
The experimental setup and procedure have been described elsewhere in detail [36,43]. In this study, the previously used experimental setup is upgraded using twolayered HeleShaw cells, with different permeability and height of the top and bottom layer. Figure 1 shows the schematic of the used setup consisting of twolayered permeability HeleShaw cells, a DIGIQUARTZ pressure transducer, a multistep programmable CHEMYX OEM syringe pump, temperature recording apparatus (three Pt100 sensors on the syringe pump, CO_{2} tank, and on top of the cell, and three EUROTHERM 2408i), CO_{2} source, imaging system, and a PC with a system controller. The cell dimensions are 50 cm × 50 cm.
Figure 1 Schematic of the experimental setup.
Two sheets of glass with dimensions of 50 cm x 40 cm and 50 cm x 10 cm are attached along their length (by glue) to produce the twolayered permeability cells. These two sheets should be placed at two different levels (representing two different gaps) for gluing. Afterward, another sheet of glass is fixed against this one, and shims of different thicknesses are placed on the sides in different layers. Small shims with the same thickness as those at the sides are placed between two glass sheets at different locations to ensure a uniform gap in different locations. The cell is tightened from four sides and center with small diagonal forces of the same magnitude.
A tank is attached to the bottom of the cell to uniformly fill the cell by water, and a stainlesssteel tank is connected to the top to provide CO_{2} (source) overlaying the water in the cell. After each test, the cell is cleaned using a mixture of distilled water and methanol. The permeability of each layer can be calculated using Equation (1),
$$k=\frac{d^2}{12}\tag{1}$$
where k and d are permeability (m^{2}) and gap (m) of each layer, respectively.
The setup is automated by maintaining a constant pressure (around atmospheric pressure) of the freephase CO_{2}, on top of the water, to prevent movement of the top boundary with more accuracy. This is achieved by checking the CO_{2} pressure at the top of the water by the pressure transducer and controlling the syringe pump based on this pressure. At each time step, when CO_{2} starts to dissolve into the water, its pressure decreases gradually. At this time, the syringe pump injects CO_{2} until reaching the initial pressure and maintains a steady pressure.
Images of the process are captured at every 10 sec, 20 sec, 1 min, and 2 min using a Canon EOS1Ds Mark II camera, connected to a PC, and a video of the convection finger movement in the cell is made. A table light is positioned behind the cell to create highquality images.
2.2 Fluid Composition and Properties
The fluids used in these experiments are a solution of 0.025 wt.% bromocresol green (SigmaAldrich, Inc.) and distilled and deionized water and CO_{2} gas at the ambient condition. The small amount of bromocresol green does not affect the properties of water. At low pH, the indicator color is yellow, and at high pH, the color is blue. Table 1 presents the thermodynamic properties of the fluids used in the experiments [36,43].
Table 1 Properties of the used fluids.
2.3 Specifications of the Twolayered Cells
Four twolayered permeability HeleShaw cells, with the specifications given in Table 2, are used for understanding the effects of vertical variation of permeability and layer thickness (in saline aquifers) on the behavior of convective mixing of CO_{2} and water. d, k, h, and R_{a} are the gap, permeability, height, and Rayleigh number of the layers, and subscripts 1 and 2 represent the top and bottom layers, respectively. The cells are oriented vertically, and their considered width and height for analysis are 0.50 m and 0.25 m, respectively (upper half of the cells). The calculations and analyses are performed for various durations before convection fingers touch the bottom boundary, and the experimental cells correspond to infinite depth aquifers with the gaswater contact situated at z = 0.
Table 2 Layered permeability HeleShaw cell models.
The behavior of convective mixing and dissolution of CO_{2} into the water in the layered permeability cells are compared with the results of the vertical homogeneous cells (with the gaps of 0.25 mm (cell H2) and 0.50 mm (cell H5) with permeabilities of 5266 Darcy and 21065 Darcy, respectively) [36].
2.4 Experimental Procedure
After filling the cell with water (with a pH of around 5.4) from the bottom, an equilibrium temperature is obtained between the environment, cell, and table light behind the cell. The inlet valve (through which CO_{2} enters the cell) and the outlet valve (through which CO_{2}/air exits from the cell) are opened simultaneously to purge the cell with CO_{2} and_{ }remove the air above the water. After few seconds (to be sure that the air inside the lines, CO_{2} tanks, and cell has been removed), the inlet and outlet valves are closed simultaneously, and the CO_{2} pump and pressure transducer valves are opened. CO_{2} is injected into the cell to maintain a constant pressure.
3. Results and Analysis
3.1 Analytical Approach
The experimental data consist of the quantitative (the amount of CO_{2} dissolved into the water) and the qualitative data (the captured images from the HeleShaw cells during the tests). The selected snapshots from the HeleShaw cell are processed, and their colors are replaced to clearly visualize the changes in the dissolved CO_{2}. In these snapshots, the blue color is pure water in equilibrium with air (pH of 5.4). The red color on top of the water shows the gaseous phase CO_{2}, and the red color in the water represents water with maximum dissolved CO_{2 }(pH of about 3.9). The other colors in between represent water with dissolved CO_{2} and acidity ranging between 3.9 to 5.4.
The amount of dissolved CO_{2} in water is calculated by calculating the number of moles of insitu gasphase CO_{2} present in the pump (first term), the tank (second term), and the cell above the water surface (third term) at each time duration, and subtracting them from the values obtained at the previous time duration to obtain the number of moles of the dissolved CO_{2} in water at each time interval:
$$(\Delta n_d)_{{i+1},i}= \frac{1}{R}\left[ \left(\frac{P_{p_{i+1}}V_{p_{i+1}}}{Z_{p_{i+1}}T_{p_{i+1}}}\frac{P_{p_i}V_{p_i}}{Z_{p_i}T_{p_i}}\right)+\left(\frac{P_{t_{i+1}}V_{t_{i+1}}}{Z_{t_{i+1}}T_{t_{i+1}}}\frac{P_{t_i}V_{t_i}}{Z_{t_i}T_{t_i}}\right)+\left(\frac{P_{c_{i+1}}V_{c_{i+1}}}{Z_{c_{i+1}}T_{c_{i+1}}}\frac{P_{c_i}V_{c_i}}{Z_{c_i}T_{c_i}}\right)\right]\tag{2}$$
In these calculations, it is assumed that:
$V_{t_{i+1}}=V_{t_i}=V_t$ , $V_{c_{i+1}}=V_{c_i}=V_c$ , $P_p=P_t=P_c=P_{{CO}_2}$ , $ T_t=T_c=T_s$ , $Z_t=Z_c=Z_{s,{CO}_2}$ , $Z_p=Z_{p,{CO}_2}$
Subscripts p, t, c, and s indicate pump, tank, cell (above the water surface in the cell), and system, respectively. In this setup and for this equation, tank (t) represents all lines, connections, valves, and stainlesssteel tank attached on top of the cell, and its volume (V_{t}) is calculated as 4.25E05 m^{3 }[36,43]. V_{c} is the cell volume above the water surface (that is full of CO_{2}) and is calculated by measuring the water level in the cell before each test. Various data (measured at each time interval of 10 s) are as follows: CO_{2} pressure in the system ($P_{{CO}_2}$), that is read by the pressure transducer and is fixed at a constant value, the pump volume (V_{p}), that represents the volume of CO_{2 }injected by the pump, the pump temperature (T_{p}) and cell temperature (T_{c}), that are read at each time interval, by Pt100 sensors attached to them. After calculating the number of moles of the dissolved CO_{2} into the water at each time interval using Equation 2, and determining the cumulative mole of dissolved CO_{2 }into the water, the calculated values are transformed from mole to kg/m^{2} unit by considering the molecular weight of CO_{2} and the crosssectional area of the contact surface. The calculated cumulative dissolved CO_{2} in kg/m^{2} (experimental value) is comparable with the diffusion equation (theoretical value) as stated in Equation 3,
$$\rm M\left(t\right)=\ 2C_0\sqrt{{Dt}/{\pi}}\tag{3}$$
Where M(t) represents the total amount of dissolved CO_{2 }per crosssectional area at time span t. The observed time of deviation of the experimental values from the theoretical values is attributed to the quantitative onset time of convection. The suitable region for the diffusion equation can be selected by examining the graph of cumulative dissolved CO_{2} vs. the square root of time, and the theoretical value (obtained from the diffusion equation) can be fitted on the experimental values (obtained from the experimental values of dissolved CO_{2}) by minimizing the root mean square differences. The quantitative onset time for convection or the time of deviation of experimental data from the theoretical data can be compared with the qualitative onset time for convection from the captured images after observing the first instabilities. The qualitative average value of the critical wavelength of convection fingers, observed after the onset of convection, is computed by selecting a centered horizontal segment (to remove the edge effect) of the CO_{2}water interface and dividing this distance by the number of convection fingers in this segment. The onset time of convection and the critical wavelength of convection fingers (calculated from the experiments conducted in the homogeneous cell) are transferred to dimensionless forms using a length scale and a time scale as given in Equations 4 and 5, and the values can be compared with the numerically predicted values [14,15,18,19,21,22,44,45,46,47]
$$L=\phi \mu D/\Delta \rho gk\tag{4}$$
$$T=L^2/D=D(\phi \mu /\Delta \rho gk)^2\tag{5}$$
Another critical parameter is the dissolution flux (after the onset of convection in the constantflux regime), which is also calculated for each test using experimental measurements.
3.2 Experimental Results and Analyses
Four tests are carried out using heterogeneous cells (twolayered cells with different heights and permeabilities), as shown in Table 2. The tests of the heterogeneous system with a 0.25 mm gap in the first layer (tests L1 and L2) are compared with a vertically homogeneous system with a uniform gap of 0.25 mm (test H2) and the tests of the heterogeneous system with 0.50 mm gap in the first layer (tests L3 and L4), are compared with a vertically homogeneous system with a uniform gap of 0.50 mm (test H5). Comprehensive results of convective mixing in the homogeneous cells are presented elsewhere [36]. The quantitative measurements of dissolved CO_{2} in water are accurately made only in the layered permeability cell L3, and in the others, only images are analyzed, and qualitative data are acquired. The videos created from the images are speeded up with the scales of 1/1600 and 1/3200 and have been uploaded to the following URL address for public awareness of the process:
https://www.youtube.com/playlist?list=PLfLgKEdPkuBPIlO1JDrmUJ0XnSCJOIBJD
3.2.1 Low Permeable Top Layer
Snapshots of convective mixing of CO_{2} and water, and the changes in the dissolved concentration of CO_{2} in water, over a certain time duration in the layered permeability cell L1 (when a low permeability layer with a gap of 0.25 mm is on top of a high permeability layer with a gap of 0.50 mm) are shown in Figure 2. This figure also compares the convective mixing process of the layered permeability cell L1 with the homogeneous cell H2 (with a gap of 0.25 mm). The first snapshots (Figure 2a) show the cells at the initial state (i.e., before the introduction of CO_{2}). After the introduction of CO_{2} on top of the cell (i.e., above the water), CO_{2} starts to dissolve into the water by diffusion, and the color of the top layers of water changes from blue (with a pH of 5.4) to yellow (with a pH of less than 5.4), because of changes in the acidity level of water by the dissolved CO_{2}. When the thickness of this diffusive layer increased sufficiently, gravitational instability occurs, and convective mixing commences at this time (onset time for convection). The calculated onset time of convection (t_{c}^{*}) in these two systems are found to be 897 sec and 970 sec using quantitative and qualitative measurements, respectively. The average size of the convection fingers (the qualitative critical wavelength of convection fingers, λ_{c}^{*}) is found to be 0.01178 m approximately. Considering the length scale (L = 9.14E05 m) and time scale (T = 4.43 sec) in the homogeneous cell H2, the dimensionless values of time corresponding to the onset of convection (t_{c}) calculated from quantitative and qualitative data are found to be 202.45 and 218.93, respectively, and the critical wavelength of convection fingers (λ_{c}) is 128.83 [36]. By comparing the snapshots of the tests carried out in the layered permeability cell and homogeneous cell at the same time, it is observed that tests of both cells display almost similar behavior before the convection finger touches the second layer (in the layered cell) at around 21591 sec (Figure 2b and 2c). This means that the top layer of the layered permeability cell is the dominant layer, which affects the behavior of convection fingers, like the onset time of convection, critical wavelength, size, and growth of convection fingers after convection starts and in the constant flux regime (before the bottom layer comes in contact with the convection fingers). Even though the convection fingers in the H2 cell proceed downward with time, the convection fingers in the L1 cell do not penetrate the high permeable bottom layer. After the convection fingers touch the interface (between two layers) at around 21591 sec in the L1 cell, they merge together to form fingers of higher wavelengths (Figure 2d, 2e, and 2f). No further downward penetration of fingers is observed.
Figure 2 Snapshots of convective mixing of CO_{2} and water in layered permeability cell L1 and homogeneous cell H2; (a) t = 0 sec, (b) t = 3391 sec, (c) t = 21591 sec, (d) t = 39853 sec, (e) t = 58034 sec, (f) t = 76274 sec.
The changes in the dissolved concentration of CO_{2} in water over a time duration in layered permeability cell L2 and homogenous cell H2 are shown in Figure 3, where a low permeability layer (with a gap of 0.25 mm) is on top of a high permeability layer (with a gap of 0.50 mm). The water height in this cell is less than the L1 cell. A comparison of the snapshots of the cell L2 with the cell H2 (at the same time) from Figure 3 shows similar results. While the behavior of L2 and H2 cells are identical initially (after starting the test), the convection fingers in the cell L2 do not enter into the bottom layer (with higher permeability and pore volume) until 76269 sec (that is the final recorded time). They grow and their wavelength increases when they touch the interface (between two layers). However, the convection fingers of the H2 cell grow and move downward.
Figure 3 Snapshots of convective mixing of CO_{2} and water in layered permeability cell L2 and homogeneous cell H2; (a) t = 0 sec, (b) t = 3405 sec, (c) t = 21669 sec, (d) t = 39788 sec, (e) t = 58029 sec, (f) t = 76269 sec.
By comparing the layered permeability cells L1 and L2, it can be stated that the behavior of convection fingers formed in the top layer (of lower permeability) when they touch the interface is not dependent on the size of the fingers (width and length). The convection finger wavelength increases without any advancement into the bottom layer because there is no sufficient mass of CO_{2} present on top of the bottom layer to continue the convective flow and the fingering process in that layer. In this condition, the diffusion mechanism dominates the convection flow. These results are consistent with previous numerical simulation results [3], which have concluded that in a layered permeability cell, the properties of the top layer (of lower permeability) only dominate the whole process till the end.
3.2.2 High Permeable Top Layer
The changes in the dissolved concentration of CO_{2} in water, over a certain time period, in a layered permeability cell L3, are displayed in Figure 4, when a high permeability layer (with a gap of 0.50 mm) is on top of a low permeability layer (with a gap of 0.25 mm). Moreover, the behavior of convection fingers in this cell and the homogeneous cell H5 (with the gap of 0.50 mm) is also compared (at the same time) in Figure 4, which shows that the behavior of both cells is almost similar before the convection fingers touch the second layer in the L3 cell, e.g., at about 771 sec and 6341 sec (Figures 4b and 4c). Similar onset time of convection, the critical wavelength, shape, position, and growth of convection fingers with time, are observed in both L3 and H5 cells. The onset time values of convection (t_{c}^{*}) in these two cells calculated using quantitative and qualitative measurements are found to be 489 sec and 230 sec, respectively, which are much lower than the calculated time values for the previous tests (where the low permeable layer was on the top). The qualitative critical wavelength of the convection fingers (the) (λ_{c}^{*}) is 0.00514 m, which is lower than the value found in the previous tests (0.01178 m). Considering the length scale (L = 2.29E05 m) and time scale (T = 0.2769 sec) in the homogeneous cell H5, the dimensionless onset time values for convection (t_{c}) calculated from the quantitative and qualitative data (t_{c}) are 1765.88 and 830.58, respectively, and the critical wavelength of convection fingers (λ_{c}) is 224.86 [36]. When the convection fingers of the L3 cell touch the bottom layer (of less permeability and pore volume), the downward movement speed of convection fingers decreases than cell H5 (Figure 4c). Then these convection fingers (formed in the top layer of higher permeability) start to merge and then penetrate the bottom layer with higher wavelengths. Although the convection fingers in the H5 cell reach the bottom boundary of the cell after 22965 sec (Figure 4f), the fingers in this layered cell can reach only till the middle of the cell at this time (Figure 4f).
Figure 4 Snapshots of convective mixing of CO_{2} and water in layered permeability cell L3 and homogeneous cell H5; (a) t = 0 sec, (b) t = 771 sec, (c) t = 6341 sec, (d) t = 11881 sec, (e) t = 17500 sec, (f) t = 22965 sec.
Figure 5 compares cumulative dissolved CO_{2} in the layered permeability cell L3 with homogeneous cells H2 and H5 with gaps of 0.25 mm and 0.50 mm, respectively. It can be seen that the rate of dissolution of CO_{2} in the L3 cell is the same as the rate of dissolution in the H5 cell with the gap of 0.50 mm until 12693 sec (≈5.986E10 kg/s), and it seems that the top layer with higher permeability is dominant until this time, while the convection fingers have touched the bottom layer at around 6341 sec. After 12693 sec, the rate of dissolution of CO_{2} into the water in the L3 cell (with k_{1}/k_{2} of 4) is reduced by 50% (to 3.152E10 kg/s) due to the effect of the bottom layer (of lower permeability). However, it is still more than twice the dissolution rate of CO_{2} in the H2 cell (1.265E10 kg/s) with the same properties of the bottom layer. It causes a considerable difference between the amount of dissolved CO_{2} into the water in L3 and H2 cells (cumulative dissolved CO_{2} of 1.9E5, 1.5E5, and 5.4E6 kg till 30000 sec in the cells H5, L3, and H2, respectively). It can be seen as that the final cumulative dissolved CO_{2} in the L3 cell is more than the H5 cell, as at the final recorded time, the convection fingers reach the bottom boundary of the H5 cell, but the convection fingers reach only till the interface (almost middle of the cell) of the L3 cell at that time. This comparison reveals the impact of the highpermeable top layer of a layered cell on increasing the rate of CO_{2} dissolution into the water and approaching the rate of CO_{2} dissolution of the layered permeability cell to the homogeneous cell H5.
Figure 5 Comparing the cumulative amount of dissolved CO_{2} in layered permeability cell L3 with homogeneous cells H2 and H5.
Figure 6 presents the changes in the dissolved concentration of CO_{2} in water over time in layered permeability cell L4, when a high permeability layer (with a gap of 0.50 mm) is on top of a low permeability layer (with a gap of 0.25 mm). The height of water in this test is lower than the layered permeability cell L3. By comparing these images with the homogeneous cell H5 (with the gap of 0.50 mm) in this figure at the same time periods, it is observed that cell L4 shows similar results to the cell L3 (mentioned in the previous section). The behaviors of layered and homogeneous cells are identical until the convection fingers touch the bottom layer at around 780 sec (Figure 6b). After the convection fingers touch the bottom layer of lower permeability and pore volume in the L4 cell, the downward movement speed of convection fingers decreases, compared to the H5 cell. The convection fingers formed in the top layer (with higher permeability) grow, and they penetrate the bottom layer with higher wavelengths. This result confirms that the behavior of convection fingers when they touch the layer interface is independent of the height of the top layer (with high permeability) and dimension (width and length) of the fingers.
Figure 6 Snapshots of convective mixing of CO_{2} and water in layered permeability cell L4 and homogeneous cell H5; (a) t = 0 sec, (b) t = 780 sec, (c) t = 6340 sec, (d) t = 11879 sec, (e) t = 17519 sec, (f) t = 22979 sec.
A previous numerical simulation study has confirmed this behavior of the convection fingers for layered permeability models with strong heterogeneity (e.g., with a permeabilityratio of 6.0). In contrast, in layered permeability models with weak heterogeneity (e.g., with a permeability ratio of 2), the top layer properties dominate the whole process [3]. While in this work, it is shown that convection fingers form fingers of higher wavelength and penetrate the bottom layer of low permeability, another simulation study has demonstrated that the convection fingers disappear in the bottom layer of low permeability, and reduce the convection effect [40].
3.3 Numerical Simulation Models
Eclipse100 flow simulator (black oil) is used to visualize the behavior of convection fingers in the layered permeability models when they reach the bottom layer [48]. A more detailed quantitative study of this problem is presented before [3]. Similar to the experimental tests, the simulation models are of twophase flow, twodimensional, and are initialized with a gas cap containing freephase CO_{2} with constant pressure on top and an aquifer with water below. This CO_{2} phase causes maximum CO_{2} concentration on top of the aquifer after the first timestep due to the dissolution of CO_{2} into the water by diffusion. Darcy and Fick's laws are considered as the governing equations in the simulation models, and the boundary conditions in simulation models are identical to those in the experiments. The thermodynamic properties of the used fluids in the experiments are given in Table 1. The size and other features of the simulation models and experimental models are the same. In the simulation models, the porosity is considered as one, and the permeabilities are changed a bit to have an equal length scale (L) and time scale (T) in the simulation and experimental models. The grid sizes are very fine, and the critical wavelength of the perturbations, which most easily gives rise to instability, can be considered as an indication of an appropriate grid size of the models in a numerical simulation. The horizontal grid block size with 1/20 of critical wavelength, obtained from linear stability analysis (by considering the higher permeability value), is a suitable size for simulation of this behavior [3,10]. The vertical grid block size equals the horizontal grid block size, and the same grid block resolution is considered in both layers. The time steps are fine enough to capture the onset time of convection with high accuracy. The simulation results in this section are based on perturbation introduced by numerical roundoff errors in the finitedifference flow simulations, and the convective mixing starts when one introduces a perturbation from the pure diffusion profile of CO_{2} into the water (below the phase contact).
Figure 7 compares the convection finger behavior in the experimental and simulated layered permeability models L1 and L3, when they encounter the bottom layer. The results of the simulated model L1 are consistent with the experimental results of cell L1, which show similar behavior of the convection fingers, where the convection fingers (formed in the top layer of low permeability) merge and form fingers of higher wavelengths after reaching the bottom layer of high permeability. In the simulated model, the fingers penetrate the bottom layer at some points with smaller wavelengths, but the penetrated fingers do not continue to grow downward and vanish. While in the experimental tests, penetration in the bottom layer by fingers is not observed. In the simulated model L3, where the top layer is of higher permeability, the finger wavelengths are seen to increase after convection fingers touch the bottom layer of lower permeability. The convection fingers are found to penetrate the bottom layer when their wavelength increases enough. Similar results are observed in the experimental layered permeability cell L3 also, which concludes that simulated and experimental results are consistent with each other.
Figure 7 Comparison of conve ction fingers in experimental and simulation layered permeability models L1 and L3.
4. Conclusions
For investigating the effect of layer properties (thickness and permeability) on the convective mixing process and dissolution of CO_{2} into the water, four tests with twolayered permeability HeleShaw cells (heterogeneous system) were carried out, and the behavior of convection fingers was analyzed and compared with the test results of homogeneous cells. The observed experimental results were consistent with the numerically simulated results. The achieved results are concluded below:
 1. In the dissolution of CO_{2} into the water in twolayered permeability geometries, the top layer controls the onset time of convection, critical wavelength of convection fingers, and early dissolution flux of CO_{2} into water.
 2. In early times, the behavior of convection fingers and the rate of CO_{2} dissolution in the layered permeability cells is similar to the homogeneous cells with the same properties of the top layer in the layered permeability cells. However, the properties of the bottom layer can affect the whole mixing process.
 3. In the case of a low permeability top layer, the convection fingers merge, and their wavelengths increase after touching the bottom layer. However, the fingers do not penetrate the bottom layer and continue the convective flow further due to insufficient mass of CO_{2} on top of the bottom layer, and the diffusion mechanism dominates the convective flow.
 4. In the case of a high permeability top layer, the convection fingers formed in the top layer form fingers of higher wavelengths after touching the bottom layer and then penetrate the bottom layer.
 5. In the layered cell, the existence of a low permeable bottom layer below a high permeable top layer reduces the dissolution rate of CO_{2} by 50% (k_{1}/k_{2} of 4) after fingers reach the bottom layer. However, this rate is still much higher (more than twice) than the dissolution rate in the homogeneous cell with the same permeability of the bottom layer. This means that the existence of a high permeable layer on top where the instability occurs is favorable for fast and early dissolution of CO_{2} into the water.
 6. The behaviors of convection fingers (when they touch the interface) are independent of the height of the layers and the size of the fingers.
Nomenclature
C_{0} 
= 
Solubility of CO_{2}, kg/m^{3} 
D 
= 
Diffusion coefficient, m^{2}/s 
d 
= 
Gap, mm 
d_{1} 
= 
Gap of top layer, mm 
d_{2} 
= 
Gap of the bottom layer, mm 
G 
= 
Horizontal gap between barriers, m 
g 
= 
Acceleration of gravity 
H 
= 
Vertical gap between barriers, m 
h 
= 
Thickness 
h_{1} 
= 
Height of top layer, m 
h_{2} 
= 
Height of bottom layer, m 
k 
= 
Permeability 
k_{1} 
= 
Permeability of top layer, D 
k_{2} 
= 
Permeability of bottom layer, D 
k_{v} 
= 
Absolute vertical permeability 
k_{s} 
= 
Barrier permeability, md 
k_{b} 
= 
Background permeability, md 
L 
= 
Length scale, m 
M(t) 
= 
Total dissolved CO_{2 }accumulated after t per crosssectional area, kg/m^{2} 
P 
= 
Pressure, bar 
R 
= 
Universal gas constant, m^{3} bar K^{−1} mol^{−1} 
S 
= 
Length of barriers, m 
T 
= 
Temperature, K; Time scale, sec 
t 
= 
Time, sec 
V 
= 
Volume, m^{3} 
Z 
= 
Gas compressibility factor 
z 
= 
Space parameter 
Δρ 
= 
Mass density increase for fully CO_{2} saturated water 
φ 
= 
Porosity 
µ 
= 
Brine viscosity 
Subscripts: 

c 
= 
Cell 
p 
= 
Pump 
s 
= 
System 
t 
= 
Tank 
Acknowledgments
This work was supported by the BIGCCS Centre under the Norwegian research program Centres for Environmentfriendly Energy Research (FME). The authors acknowledge the following partners for their contributions: Gassco, Shell, Statoil, TOTAL, GDF SUEZ and the Research Council of Norway (193816/S60).
Author Contributions
Conceptualization: Amir Taheri, Erik Lindeberg, Dag WesselBerg; Methodology: Amir Taheri, Erik Lindeberg, Nanji J. Hadia; Validation: Amir Taheri, Erik Lindeberg, Dag Wessel Berg; Analysis: Amir Taheri, Erik Lindeberg, Dag Wessel Berg; Investigation: Amir Taheri; Supervision: Ole Torsæter, Dag WesselBerg; Project adminisration: Ole Torsæter, Dag WesselBerg; Funding acquisition: Ole Torsæter, Dag WesselBerg.
Competing Interests
The authors have declared that no competing interests exist.
References
 Lindeberg E, WesselBerg D. Vertical convection in an aquifer column under a gas cap of CO_{2}. Energy Convers manag. 1997; 38: S229S234. [CrossRef]
 EnnisKing J, Paterson L. Role of convective mixing in the longterm storage of carbon dioxide in deep saline formations. Proceeding of SPE Annual Technical Conference and Exhibition; 2003 October 58; Denver, CO, USA. Richardson, TX: OnePetro. [CrossRef]
 Taheri A, WesselBerg D, Torsaeter O, Soroush M. The effects of anisotropy and heterogeneity on CO_{2} dissolution in deep saline aquifers. Proceeding of Carbon Management Technology Conference; 2012 February 79; Orlando, FL, USA. Richardson, TX: OnePetro. [CrossRef]
 Lindeberg E, Bergmo P. The longterm fate of CO_{2} injected into an aquifer. Greenh Gas Control Technol. 2003; 1: 489494. [CrossRef]
 Lindeberg E, Vuillaume JF, Ghaderi A. Determination of the CO_{2} storage capacity of the Utsira formation. Energy Procedia. 2009; 1: 27772784. [CrossRef]
 Farajzadeh R, Salimi H, Zitha PL, Bruining H. Numerical simulation of densitydriven natural convection in porous media with application for CO_{2} injection projects. Int J Heat Mass Transf. 2007; 50: 50545064. [CrossRef]
 Pruess K, Zhang K. Numerical modeling studies of the dissolutiondiffusionconvection processduring CO_{2} storage in saline aquifers. Berkeley, CA: Lawrence Berkeley National Lab.(LBNL); 2008; LBNL1243E. [CrossRef]
 Green CP, EnnisKing J. Effect of vertical heterogeneity on longterm migration of CO_{2} in saline formations. Transp Porous Media. 2010; 82: 3147. [CrossRef]
 Pau GS, Bell JB, Pruess K, Almgren AS, Lijewski MJ, Zhang K. Highresolution simulation and characterization of densitydriven flow in CO_{2} storage in saline aquifers. Adv Water Resour. 2010; 33: 443455. [CrossRef]
 Lindeberg E, WesselBerg D. Upscaling studies of diffusion induced convection in homogeneous and heterogeneous aquifers. Energy Procedia. 2011; 4: 39273934. [CrossRef]
 Farajzadeh R, Ranganathan P, Zitha PL, Bruining J. The effect of heterogeneity on the character of densitydriven natural convection of CO_{2} overlying a brine layer. Adv Water Resour. 2011; 34: 327339. [CrossRef]
 Vosper H, Kirk K, Rochelle C, Noy D, Chadwick A. Does numerical modelling of the onset of dissolutionconvection reliably reproduce this key stabilization process in CO_{2} Storage? Energy Procedia. 2014; 63: 53415348. [CrossRef]
 Caltagirone JP. Stability of a saturated porous layer subject to a sudden rise in surface temperature: Comparison between the linear and energy methods. Q J Mech Appl Math. 1980; 33: 4758. [CrossRef]
 EnnisKing J, Preston I, Paterson L. Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Physi Fluids. 2005; 17: 084107. [CrossRef]
 Yoon DY, Choi CK. Thermal convection in a saturated porous medium subjected to isothermal heating. Korean J Chem Eng. 1989; 6: 144149. [CrossRef]
 Tan KK, Thorpe RB. The onset of convection caused by buoyancy during transient heat conduction in deep fluids. Chem Eng Scie. 1996; 51: 41274136. [CrossRef]
 Tan KK, Sam T, Jamaludin H. The onset of transient convection in bottom heated porous media. Inte J Heat Mass Transf. 2003; 46: 28572873. [CrossRef]
 Xu X, Chen S, Zhang D. Convective stability analysis of the longterm storage of carbon dioxide in deep saline aquifers. Adv Water Resour. 2006; 29: 397407. [CrossRef]
 Riaz A, Hesse M, Tchelepi H, Orr F. Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. J Fluid Mech. 2006; 548: 87111. [CrossRef]
 WesselBerg D. On a Linear stability problem related to underground CO_{2} storage. SIAM J Appl Math. 2009; 70: 12191238. [CrossRef]
 WesselBerg D. The gravitational instability of a diffusive boundary layer; towards a theoretical minimum for time of onset of convecti. Proceeding of ECMOR XIII13th European Conference on the Mathematics of Oil Recovery; 2012 September 1013; Biarritz, France. Houten, The Netherlands: European Association of Geoscientists & Engineers. [CrossRef]
 Taheri A, WesselBerg D, Torsæter O. Prediction of the minimum onset time for convection in a gravitationally unstable diffusive boundary layer using a finite difference pressure solver. Proceeding of EUROPEC 2015; 2015 June 14; Madrid, Spain. Richardson, TX: OnePetro. [CrossRef]
 Okhotsimskii A, Hozawa M. Schlieren visualization of natural convection in binary gasliquid systems. Chem Eng Sci. 1998; 53: 25472573. [CrossRef]
 Arendt B, Dittmar D, Eggers R. Interaction of interfacial convection and mass transfer effects in the system CO_{2}water. Int J Heat Mass Transf. 2004; 47: 36493657. [CrossRef]
 Yang C, Gu Y. Accelerated mass transfer of CO_{2} in reservoir brine due to densitydriven natural convection at high pressures and elevated temperatures. Ind Eng Chem Res. 2006; 45: 24302436. [CrossRef]
 Farajzadeh R, Barati A, Delil HA, Bruining J, Zitha PL. Mass transfer of CO_{2} into water and surfactant solutions. Pet Sci Technol. 2007; 25: 14931511. [CrossRef]
 Farajzadeh R, Delil HA, Zitha PLJ, Bruining J. Enhanced mass transfer of CO_{2} into water and oil by natural convection. Proceeding of EUROPEC/EAGE Conference and Exhibition; 2007 June 1114; London, UK. Richardson, TX: OnePetro. [CrossRef]
 Farajzadeh R, Zitha PL, Bruining J. Enhanced mass transfer of CO_{2} into water: Experiment and modeling. Proceeding of EUROPEC/EAGE Conference and Exhibition; 2009 June 0811; Amsterdam, The Netherlands. Richardson, TX: OnePetro. [CrossRef]
 Neufeld JA, Hesse MA, Riaz A, Hallworth MA, Tchelepi HA, Huppert HE. Convective dissolution of carbon dioxide in saline aquifers. Geophys Res Lett. 2010; 37: L22404. [CrossRef]
 Kneafsey TJ, Pruess K. Laboratory flow experiments for visualizing carbon dioxideinduced, densitydriven brine convection. Transp Porous Media. 2010; 82: 123139. [CrossRef]
 Kneafsey TJ, Pruess K. Laboratory experiments and numerical simulation studies of convectively enhanced carbon dioxide dissolution. Energy Procedia. 2011; 4: 51145121. [CrossRef]
 Kilpatrick A, Rochelle C, Noy D. Experimental visualisation and modelling of the formation and migration of density plumes during CO_{2} storage. Flows and mechanics in natural porous media from pore to field scale. Pore2Field. 2011.
 Khosrokhavar R, Elsinga G, Mojaddam A, Farajzadeh R, Bruining J. Visualization of natural convection flow of super critical CO_{2} in water by applying schlieren method. Proceeding of SPE EUROPEC/EAGE Annual Conference and Exhibition; 2011 May 2326; Vienna, Austria. Richardson, TX: OnePetro. [CrossRef]
 Taheri A, Torsaeter O, WesselBerg D, Soroush M. Experimental and simulation studies of densitydrivenconvection mixing in a HeleShaw geometry with application for CO_{2} sequestration in brine aquifers. Proceeding of SPE Europec/EAGE Annual Conference; 2012 June 47; Copenhagen, Denmark. Richardson, TX: OnePetro. [CrossRef]
 Faisal TF, Chevalier S, Sassi M. Experimental and numerical studies of density driven natural convection in saturated porous media with application to CO_{2} geological storage. Energy Procedia. 2013; 37: 53235330. [CrossRef]
 Taheri A, Lindeberg E, Torsæter O, WesselBerg D. Qualitative and quantitative experimental study of convective mixing process during storage of CO_{2} in homogeneous saline aquifers. Int J Greenh Gas Control. 2017; 66: 159176. [CrossRef]
 Nield DA, Bejan A. Convection in porous media. Cham: Springer; 2017. [CrossRef]
 Mal'kovskii V, Pek A. Onset conditions of free thermal convection of a singlephase fluid in a horizontal porous layer with depthdependent permeability. Izv Phys Solid Earth. 1999; 35: 990994.
 Green CP, EnnisKing J. Steady dissolution rate due to convective mixing in anisotropic porous media. Adv Water Resour. 2014; 73: 6573. [CrossRef]
 Farajzadeh R, Zinati FF, Zitha P, Bruining J. Densitydriven natural convection in dual layered and anisotropic porous media with application for CO_{2} injection project. Proceeding of ECMOR XI11th European Conference on the Mathematics of Oil Recovery; 2008 September 0811; Bergen, Norway. Houten, The Netherlands: European Association of Geoscientists & Engineers. [CrossRef]
 McKIBBIN R, Tyvand PA. Anisotropic modelling of thermal convection in multilayered porous media. J Fluid Mech. 1982; 118: 315339. [CrossRef]
 Post V, Simmons CT. Free convective controls on sequestration of salts into lowpermeability strata: insights from sand tank laboratory experiments and numerical modelling. Hydrogeol J. 2010; 18: 3954. [CrossRef]
 Aggelopoulos CA, Tsakiroglou CD. Effects of microheterogeneity and hydrodynamic dispersion on the dissolution rate of carbon dioxide in watersaturated porous media. Int J Greenhouse Gas Control. 2012; 10: 341350. [CrossRef]
 Agartan E, Trevisan L, Cihan A, Birkholzer J, Zhou Q, Illangasekare TH. Experimental study on effects of geologic heterogeneity in enhancing dissolution trapping of supercritical CO_{2}. Water Resour Res. 2015; 51: 16351648. [CrossRef]
 Taheri A, Torsæter O, Lindeberg E, Hadia NJ, WesselBerg D. Qualitative and quantitative experimental study of convective mixing process during storage of CO2 in heterogeneous saline aquifers. Int J Greenh Gas Control. 2018; 71: 212226. [CrossRef]
 Kim SY, Kuznetsov AV. Optimization of pinfin heat sinks using anisotropic local thermal nonequilibrium porous model in a jet impinging channel. Numer Heat Transf A. 2003; 44: 771787. [CrossRef]
 Hassanzadeh H, PooladiDarvish M, Keith DW. Stability of a fluid in a horizontal saturated porous layer: Effect of nonlinear concentration profile, initial, and boundary conditions. Transp Porous Media. 2006; 65: 193211. [CrossRef]
 Selim A, Rees DAS. The stability of a developing thermal front in a porous medium. I Linear theory. J Porous Media. 2007; 10: 116. [CrossRef]
 Vadász P. Emerging topics in heat and mass transfer in porous media: From bioengineering and microelectronics to nanotechnology. Dordrecht: Springer; 2008. [CrossRef]
 Schlumberger. Eclipse 100 Reservoir Simulator. V 2014.1. 2014. Available from: http://www.ipt.ntnu.no/~kleppe/TPG4150/EclipseReferenceManual.pdf.