Effects of Permeability Variations on CO2 Convection in Anisotropic and Heterogeneous Saline Formations

This study simulated the natural convection of dissolved carbon dioxide (CO2) in a small-scale heterogeneous saline formation using the state module ECO2N equation in the TOUGHREACT model. A one-way downscaling approach that involves using a series of sub-models in simulation procedures was proposed to efficiently simulate problems with high-scale discrepancies. This study evaluated the effects of different degrees of small-scale permeability variations on the vertical migration of dissolved CO2. The sequential Gaussian simulation model was used to generate unconditional random permeability fields for different natural logarithm of permeability (lnk) variations (i.e., lnk variances and correlations in x and z directions). The results showed an identical transition zone of dissolved CO2 near the top boundary, where a constant CO2 gas saturation was specified. The local permeability variations can trigger fingerings and enhance the vertical convection of the dissolved CO2. The number of fingerings depends on the variations in permeability near the front interface of the dissolved CO2 (i.e., the bottom edge of the transition zone for the dissolved CO2). However, the fingering patterns and developments are constrained by the permeability variations along the fingering paths. At the same mean lnk permeability the convection fluxes increase with an increase in lnk variances. However, an increase in lateral correlations (i.e., increase in the correlation lengths in the x direction) can slightly reduce the convection fluxes at the same lnk variance. The highly variable flux rates of the dissolved CO2 occur early and the variations in the flux rate decrease with time.


IntrOduCtIOn
Carbon dioxide (CO 2 ) geological sequestration in deep saline formations is a feasible approach for mitigating CO 2 emissions. This approach involves the direct injection of captured CO 2 into compatible deep saline formations for long-term storage (Bachu and Adams 2003;IPCC 2005;Michael et al. 2010;Rosenbauer and Thomas 2010;Chasset et al. 2011;Hatzignatiou et al. 2011;Kolenković et al. 2013;Ran and Zohar 2013;Trémosa et al. 2014). The density of supercritical CO 2 is lower than that of ambient saline water and this can lead to the migration of a CO 2 gaseous phase plume upward because of the buoyancy-driven force. The injected CO 2 dissolves later in the saline water present in the saline formations and the soluble CO 2 molecules diffuse within the saline water. CO 2 dissolution increases the den-sity of saline water and triggers vertical convection in the saline formation. The speed of this density-driven process is higher than that of the dissolved CO 2 plume diffusion. Such density-driven convection is the main mechanism for controlling the migration of an aqueous CO 2 plume and enhancing dissolution in the vertical direction. Researchers have indicated that density-driven convection plays a crucial role in reducing the risk of CO 2 leakage into the surface (Lindeberg and Wessel-Berg 1997;Hassanzadeh et al. 2005;Xu et al. 2006;Pruess and Zhang 2008;Pau et al. 2010;Farajzadeh et al. 2011;Ranganathan et al. 2012).
Several studies have focused on CO 2 density-driven convection in homogeneous porous media and its relevant applications (Simmons and Narayan 1997;Hassanzadeh et al. 2005Hassanzadeh et al. , 2007Farajzadeh et al. 2007;Lu and Lichtner 2007;Kandaswamy and Eswaramurthi 2008;Pau et al. 2010; Musuuza et al. 2011). Typical approaches for creating instability convection rely on applying initial perturbations to trigger the developments of fingerings for dissolved CO 2 . Hassanzadeh et al. (2005) modeled convective mixing during the onset of instability in a homogeneous formation for CO 2 storage. Their study revealed that placing a high-density fluid above a low-density fluid might forcedly lead to the convection mechanism. Moreover, their study clarified that the CO 2 dissolution rate increases according to the densitydriven flow. Hassanzadeh et al. (2007) performed a linear stability analysis with scaling behavior to examine the onset of natural convection in a homogeneous conceptual frame. Their investigation showed that the scaling relationship effect can be considered a natural convection phenomenon. Pau et al. (2010) developed a numerical model for analyzing the instability of density-driven convection during CO 2 sequestration in homogeneous media. The results from a total velocity analysis showed that long-term CO 2 storage behavior is highly dependent on the mean permeability and density-driven convection.
The hypothetical assumption of homogeneity in formations of interest is the basis for understanding the convection dynamics of dissolved CO 2 . However, natural geological formations typically involve different degrees of heterogeneity in material permeability, and such variations can considerably influence density-driven natural convection (Royer and Flores 1994;Flett et al. 2007;Lindeberg and Wessel-Berg 2011). Natural convection influenced by permeability variations has attracted considerable attention (e.g., Riaz et al. 2006;Green et al. 2009;Farajzadeh et al. 2011;Ranganathan et al. 2012). Riaz et al. (2006) indicated that the critical time for convection fingers can occur from 2000 years to approximately 10 days for permeability variations of 1 -3000 mD. Green et al. (2009) studied the vertical heterogeneity effect on longterm CO 2 migration in saline formations. In comparison to heterogeneous shale distributions, a homogeneous medium with equivalent effective vertical permeability leads to a longer breakthrough time for deep injection and a longer onset time for convection. Farajzadeh et al. (2011) focused on the permeability heterogeneity effects on density-driven natural convection. Stochastic random permeability fields were generated based on the Dykstra-Parson coefficient to simulate natural convection and estimate the CO 2 dissolution rate in brine. Using different degrees of permeability variations and spatial correlation lengths, they found that the CO 2 dissolution rate was higher and that dissolution occurred sooner in heterogeneous media than in homogeneous media. The permeability field structure strongly dominates convection. Ranganathan et al. (2012) conducted numerical simulations to assess natural CO 2 convective processes in heterogeneous media. Different simulation conditions such as domain sizes, boundary conditions and random permeability field distributions were applied in their study. The results indicated that an increase in the permeability field variance leads to an increase and decrease in instabilities at earlier and later times, respectively. Dissolution, diffusion and convection during CO 2 storage in saline aquifers is a complex mechanism for numerical simulations. The modeling scale indirectly affects the convection behavior of dissolved CO 2 . Particularly, the number of large-scale dominant fingers decreases as the domain height increases and the increase in domain width reduces the finger flow interactions. Zhang et al. (2011) investigated the convection process using a two-dimensional model based on the mineralogical composition in the Songliao Basin, China. They used the hydrogeological parameters used by Zhang et al. (2009). Numerous cases, including an increase in the vertical permeability vector, were conducted to assess convection processes in anisotropic media. The results showed that an increase in the vertical permeability vector more strongly affects the convection process than does an increase in horizontal permeability.
CO 2 convection investigations often use relatively small simulation domains because high resolution is required for resolving small-scale fingers or permeability variations in a relatively short duration. For a large-scale problem, the simulation domain might reach hundreds of meters to several kilometers and the temporal resolution can be from several years to tens of years. This is technically inefficient for CO 2 convection numerical simulations for such a large-scale problem. We propose a one-way downscaling approximation for integrating large-scale CO 2 migration and smallscale CO 2 convection using the TOUGHREACT model and the associated ECO2N module. This downscaling approximation is expected to physically connect small-scale CO 2 natural convection simulations to the large-scale model. Figure 1 shows the general one-way downscaling approximation concept for this study. Similar to the concept used in Ni et al. (2010) and other studies Hill 2002, 2004;Afshari et al. 2008), we developed a data extraction program for collecting the simulation results from a largescale (parent) model for specified areas (small-scale or submodels) and times. The collected data include solutions for all flow and transport to the large-scale model. The solutions along the boundaries of the specified areas (i.e., the smallscale models) can be reformulated to desired values and types of boundary conditions depending on the simulation requirements. Each small-scale model can then be simulated individually using the same simulator [i.e., the TOUGHREACT/ ECO2N model; (Pruess 2005;Xu et al. 2008)].

COnCEPtuAl mOdElS And numErICAl ExPErImEntS
We illustrated the concept and natural convection process using a two-dimensional profile model similar to that used in Kongsjorden et al. (1997). This study focused on exploring the small-scale physical phenomena involved in dissolved CO 2 transport such as fingering, channeling and lateral mixing for anisotropic and heterogeneous saline formations. Although the TOUGHREACT model was developed to consider complex chemical reactions, we simply considered large-scale CO 2 migration and small-scale CO 2 convection for illustration purposes. To present numerical examples for small-scale permeability variations this study used the sequential Gaussian simulation (SGSIM) algorithm (Deutsch and Journel 1998) to generate unconditional random permeability fields with various variances and correlation lengths. The SGSIM permeability fields were then integrated with the numerical simulator TOUGHREACT/ECO2N to systematically model the large-scale migration and small-scale convection of CO 2 for the offshore gas field Sleipner.
Since September 1996 the offshore gas field Sleipner, located in the middle North Sea, has injected 1 Mt CO 2 year -1 (Kongsjorden et al. 1997). These field operations involve injecting CO 2 into a saline formation called the Utsira formation, which is a sandstone formation that consists of homogeneous, high-permeability, fine-grained sand with microfossil fragments deposited on a shallow marine shelf. The sand formation is overlaid by thick Hordaland shale. The top of the formation is located at approximately 800 m below the sea bottom and the formation thickness varies from 150 -250 m. The pressure in the formation is nearly hydrostatic and the temperature is estimated to be 25 -28°C (Korbøl and Kaddour 1995).

large-Scale models
The gas field Sleipner conceptual model has been used in previous CO 2 simulation investigations (Torp and Gale 2004;Cavanagh 2013;Chadwick and Eiken 2013;Deflandre et al. 2013;Fornel and Estublier 2013;Cavanagh and Haszeldine 2014). We used the conceptual model proposed by Pruess (2005) for the two-dimensional large-scale model. Figure 2 shows the modified conceptual models used for CO 2 migration and convection simulations. The 6000 m × 180 m simulation domain contains five reservoirs composed of relatively high-permeability sandstone of 3 × 10 -12 m 2 . The simplified conceptual model consists of four thin, low-permeability layers composed of shale stone. These low-permeability layers are considered cap rocks acting as barriers preventing upward CO 2 migration. The permeability value for the cap rocks was fixed at 1 × 10 -14 m 2 . No heat or mass was allowed to pass through the boundaries. The top and bottom flow simulation boundaries were assigned as impermeable boundaries. The right boundary was fixed at a specific hydrostatic pressure. The left boundary was considered a symmetry plane and a no flow boundary condition was specified, except for the injection location where the pressure was known before the injection. To prepare the initial condition for the large-scale model simulation, pressure is held constant at 1.1 × 10 7 kPa at the injection node elevation (22 m) and the system was run to gravity equilibrium.
The injection point was located in the lowest layer of the simulation domain (Fig. 2). The CO 2 injection rate was fixed at 0.1585 kg s -1 . This condition represents a horizontal injection well. A 2-year injection was applied to an area of 1 m × 1 m. According to the hydrostatic pressure applied in the simulation domain for initial conditions, the total simulation time was 62 years, including the 2-year injection and 60-year migration after the injection is stopped. Table 1 summarizes the boundary conditions and the associated hydrogeological parameters for the simulations of flow and CO 2 injection.

downscaling models
The one-way downscaling approach involves using a series of interpolations for extracting information from the large-scale model. Figure 3 shows the conceptual framework of the downscaling approach. A similar concept has been widely used in many studies on problems with high-scale discrepancies Hill 2002, 2004;Afshari et al. 2008).
The hierarchy tree begins from a main model and fol-lows a series of sub-models. Different sub-model layers bridge the relationships between the parent and child models. A parent model can have many sub-models; however, a child model has only one parent model because the information source is identical. We defined the sub-models extracted from the main model as first-level sub-models, Relative permeability (liquid) van Genuchten which were denoted as "sub 1" or "sub 2." The sub-model extracted from sub-model 1 is the second-level sub-model, which is denoted as "sub 1-1" or "sub 1-2." The same rule can be applied to other sub-models in the hierarchy tree. Figure 4 shows the downscaling processes corresponding to the modeling structure shown in Fig. 3. The smallest submodel (i.e., sub 1-1-1) in Fig. 4 shows the result obtained by interpolating solutions from the main model and two intermediate sub-models (i.e., sub 1 and 1-1). Any sub-model can be a new main model for conducting a numerical simulation based on the information from its parent model. The information transferred from the large-scale model to the small-scale models might introduce interpolation errors in connecting different sub-model levels. The number of levels and scale discrepancies between sub-models must be adjusted depending on the problems of interest. For practical problems the scale of a sub-model must cover at least three to five numerical cells in the parent model to accurately capture the hydrodynamic variations in the parent model. Novel features were developed for the downscaling model in this study. The interpolations in the downscaling model are performed to extract solutions (i.e., pressure, CO 2 gas saturation SG, or XCO2a) and parameters (i.e., permeability, porosity, and others) from the main model and a series of sub-models. With predefined sub-model sizes and numbers of cells in different directions, the developed search loops can extract solutions from the TOUGHREACT model input and output. This study simply used the bilinear interpolation algorithm to obtain the sub-model information. Other interpolation algorithms can be implemented using the same modeling procedures. The downscaling model then integrates the extracted solutions and parameters to become the input files for new individual sub-models. This feature enabled us to conduct simulations using a single TOUGHREACT model. The sub-model outputs for specified times were determined according to the output commands used in the main model.

numerical Examples of dissolved CO 2 Convection
We focused on the convection behaviors of dissolved CO 2 influenced by small-scale permeability variations. To evaluate CO 2 convection caused by permeability variation the simulation domain was substantially reduced to handle the small-scale permeability variations and resolve the detailed fingerings involved in convection processes. In numerical examples we used a modeling size of 5 m × 5 m, which is similar to that used in Hassanzadeh et al. (2005) and Pruess and Zhang (2008). At the simulation scale the modeling size must be three orders of magnitude smaller than that of the original main model shown in Fig. 3. Following the procedures described in the previous section we used the downscaling approach to obtain a submodel that fitted the scale of our numerical examples (i.e., 5 m × 5 m models). The sub-model (i.e., sub 1-1-1) in Fig. 4 was the basis for conducting numerical assessments in this study. Figure 5 shows the conceptual model and a random permeability field for the numerical examples. The origin of the coordinate was adjusted to (0, 0) for to enhance the presentation. To evaluate the permeability heterogeneity effect on natural dissolved CO 2 convection we maintained only the boundary SG values along the top boundary for the selected sub-model and assumed the initial XCO2a and SG concentration to be zero for the entire sub-model. However, the initial condition for flow in the sub-model was fixed to the interpolated pressure from the parent model (i.e., sub 1-1). The hydrostatic boundary condition was specified on the left and right boundaries. Table 1 lists the boundary conditions and the associated hydrogeological parameters for the conceptual model for the test cases. All modeling cases were maintained under the same hydrogeological conditions except for the variation in the small-scale permeability. The temperature gradient was not considered in the simulations and was fixed at 37°C throughout the entire modeling domain in each case. In addition the salinity was fixed at 3.2% in numerical cases to evaluate solely the permeability variation effect on dissolved CO 2 migration. This is different from previous investigations in which the variations in salinity could directly influence CO 2 solubility. High salinity can reduce CO 2 solubility and reduce the convection driving force (Zhang et al. 2011).
This study employed the Gaussian-based random field (RF) generator, an SGSIM technique in Geostatistical Library, to generate permeability fields for the test cases. The exponential covariance model for the natural logarithm of permeability (lnk) was specified to generate the RFs. Table 2 summarizes the geostatistical RF parameters for different cases. The lnk variances ( lnk 2 v ) and anisotropic ratios ( r m ) represent the degrees of spatial lnk variations and directional correlations, respectively. Case H in Table 2 shows the homogeneous case used for comparison purposes. However, Case 2 was the base case for comparisons of heterogeneous cases. We were interested in the XCO2a fingerings triggered by different degrees of permeability variations. In addition, the quantifications of the fluxes caused by different cases are systematically presented in the following section.

rESultS And dISCuSSIOn
3.1 CO 2 migration in large Scale model Figure 6 shows the simulation results for the largescale model at 2, 4, 22, and 62 years. The CO 2 gas saturation (SG) distribution showed that the injected CO 2 migrated to the second layer (left column in Fig. 6). The aqueous phase CO 2 distribution and a dissolved CO 2 mass fraction (XCO2a) showed the same patterns. However, a relatively high mass aqueous phase CO 2 fraction was obtained under the conditions used in this study (right column in Fig. 6). Because of the low-permeability cap rocks, most gaseous phase CO 2 was trapped in the injected layer during the injection. Parts of the gaseous phase CO 2 passed through the cap rocks and were concentrated under the cap rocks in the upper layers. After the injection was stopped no pressure stress source was provided in the system. The gaseous phase CO 2 plumes under the cap rocks gradually decreased with time and moved laterally along the bottom of the cap rocks. Figure 7 shows the temporal variations in SG and XCO2a along the bottom of the cap rock located above the injection layer. The spatial variations in SG (Fig. 7) clearly showed systematical trends based on the SG concentration   front after 2 years. The SG front after 2 years reached approximately 2000 m from the injection point. After the injection was stopped the SG between the concentration front and the injection point decreased with time. The SG distributions increased with time at distances greater than 2000 m (i.e., the concentration front) when the SG distributions were compared with that after 2 years. Therefore, such variations in the SG distribution at distances > 200 m lead to a piston-like movement for the XCO2a front (Fig. 7b). The density-driven behavior after CO 2 injection has been addressed and identified as soluble trapping in many previous studies (Hassan-zadeh et al. 2005;Xu et al. 2006Xu et al. , 2008Pruess and Zhang 2008;Pau et al. 2010;Ranganathan et al. 2012). Figures 8 and 9 show the dissolved CO 2 mass fraction distributions in Cases H and 1, respectively. The results clearly show similar transition zones from the upper boundary to approximately z = 4.5 m. The transition zone result was applied for all other cases in this study. We found that the transition zone was significantly influenced by salinity in the modeling domain. Regardless of the mesh size, pressure distribution and SG value along the top boundary, the transition zone was the same (i.e., approximately down to the z = 4.5 m profile). Figure 8 shows a comparison of two mesh sizes for the dissolved CO 2 mass fraction simulations in the homogeneous case. Figure 8a shows the results for a mesh size of 10 cm × 10 cm, and Fig. 8b shows the results for a mesh size of 6.25 cm × 6.25 cm. The differences in dissolved CO 2 migration in Fig. 8 reveal that the XCO2a fingerings can be triggered by numerical instability. In Fig. 8a relatively uniform fingerings were obtained after 200 days. This result guided us to consider the 6.25 cm × 6.25 cm mesh size to avoid the numerical instability effect on evaluating CO 2 convection. The 6.25 cm × 6.25 cm mesh sixe was then used for other cases in this study. At the selected mesh size, flow and transport instability in cases other than Case H can be considered to be solely influenced by the permeability variations. The fingering initiation and development mechanisms differed between homogeneous (Fig. 8a) and heterogeneous ( Fig. 9) cases. Although a relatively small variation ( lnk 2 v = 0.5) was applied to the random permeability field, the fingerings in the heterogeneous case (Case 1) began developing after 10 days. The fingering fronts in the heterogeneous case (Case 1) touched the bottom boundary after 300 days. In Fig. 9 the irregularly distributed fingers show patterns and sizes different from those induced by numerical instability. In the early time periods the fingering fronts moved down to the center of the modeling domain in less than 50 days. In addition to the vertical movement the fingerings began to expand laterally after 100 days. The continuously lateral expansions merged fingers after 200 days (Fig. 9). The results are in agreement with those reported by previous investigations, indicating that formation heterogeneity can enhance dissolved CO 2 convection (Green et al. 2009;Farajzadeh et al. 2011;Ranganathan et al. 2012).

lnk Variance and Anisotropic ratio Effects on CO 2 Convection
We further selected seven cases to assess the effects of spatial lnk variability on CO 2 convection. Case 2 was used as the base case for assessing the convection behaviors under small-scale formation heterogeneity. The lnk 2 v values were specified as 0.5, 1.5, and 2.0 in Cases 1, 3, and 4, respectively, and the correlation lengths in the x and z directions (i.e., x m and z m , respectively) for these cases were fixed at 1.0 m. In Cases 5, 6, and 7, x m was changed to produce different anisotropic ratios (  dissolved CO 2 mass flux along the z = 4.5 m profile, which is immediately below the interface where the fingerings begin to develop. Again, all cases had similar transition zones near the top boundaries. The results in Fig. 10 clearly show similar patterns because the random seed for generating the RF was the same for each heterogeneous case. The increase in lnk 2 v can lead to an increase in the spread of XCO2a distributions. The fingering in highlnk 2 v cases, such as Cases 3 and 4, tended to spread and merge early with nearby fingerings (see the fingerings between x = 4 m and x = 5 m after 30 days). In Cases 1 and 2, two fingerings merged after 100 days. Similar behaviors were observed for two fingerings between x = 0 m and x = 1 m. The two fingerings merged after 100 days in cases with lnk 2 v > 1.0. In addition a high lnk 2 v value enhances the movement of fingering fronts. However, the differences were not significant in the early time periods.
We noted that the SG distributions in transition zones varied from 0.4% (top boundaries) to 0.04% (the interfaces of fingerings). In this transition zone, XCO2a was not influenced by the small-scale lnk variations. Below the lower boundary of the XCO2a transition zone (approximately along the z = 4.5 m profile), the small-scale lnk variation effects on XCO2a distributions were significant depending on the permeability values near the interfaces (see Figs. 11a and b). The coexistence of local high lnk values and high mass fluxes revealed a direct correlation between local high lnk and fingering initiations. Furthermore, we calculated the mass flux rates of the dissolved CO 2 at specified times. Figures 11b to f indicates the dissolved CO 2 mass flux rate for the specified times along the A -A' profile (z = 4.5 m). The simulation results showed that the dissolved CO 2 mass flux rates along the z = 4.5 m profile were highly variable in the beginning. However, the dissolved CO 2 mass flux rates decreased with time. Fingering occurs according to the locations of relatively high permeability values. Only relatively high permeability values near the interfaces of the fingerings triggered initial fingerings at the early time points (Figs. 11a and b). However, the development and mass flux for a fingering might not be directly correlated with the high permeability values near the lower boundary of the XCO2a transition zone. The permeability below the fingerings and along the fingering paths substantially controlled the devel-opment and the dissolved CO 2 mass flux of a fingering at different times. Figures 12 and 13 show results comparison based on different lnk anisotropic ratios. In the comparison a fixed lnk 2 v of 1.0 and y-direction correlation length of 0.3 m were used. In Cases 5, 6, and 7 the x-direction correlation lengths were increased to 0.9, 1.5, and 3.0 m, respectively. The results in Fig. 12 show that an increase in the x-direction correlation lengths can reduce fingering development in both lateral spreading and front movements (Fig. 12a). After 200 days the low anisotropic ratios caused XCO2a values to be high in the center of fingerings near the lower boundary of the XCO2a transition zone (Fig. 12b). However, the fingering sizes for high anisotropic ratios were obviously smaller than those in the cases with low anisotropic ratios. Figure 13 shows a comparison of the dissolved CO 2 mass fluxes for different anisotropic ratios at different times. The permeability k profile in Fig. 13a clearly shows relatively small k variations in cases with large anisotropic ratios. Such small k variations along a profile can result in small CO 2 convections. The permeability and dissolved CO 2 mass fluxes along the A -A' profile showed consistent behaviors similar to those in Fig. 11. Fingering occurrences were based on the relatively high permeability near the lower boundary of the XCO2a transition zone (i.e., approximately along the z = 4.5 m profile). The dissolved CO 2 mass flux decreased with time. At early time periods a decrease in the anisotropic ratio reduced the dissolved CO 2 mass flux. In addition, the differences in the anisotropic effects on mass flux became minor at late time points.

Effects of mean Permeability
Cases 8 and 9 were used to investigate the mean permeability variation effects on the migration of dissolved CO 2 . Table 2 lists the relevant hydrogeological parameters for these two cases. Figure 14 shows a comparison of the XCO2a distributions in Cases 8 and 9 at 30 and 200 days, respectively. When the statistical structure of the RFs was applied to Cases 8 and 9 (Table 2), the fingering development was significant in Case 8 compared with the fingerings in Cases 2 and 9. In Case 8 the fingerings reached the bottom boundary after 30 days. However, in Case 9, the initial fingerings began to develop after 200 days. Even when fingerings were triggered by local permeability variations, low mean permeability reduced the CO 2 convection activity. We conclude that mean permeability dominates the overall CO 2 convection in reservoir formation. However, the initiation, development and lateral mixing of fingerings were controlled significantly by the permeability variations near the XCO2a transition zone lower boundary and along the fingering paths. Again, the transition zones were similar because the salinity was constant in the simulation domain.

Comparisons of mass Fluxes Among Cases
The mass convection rates were calculated in this study based on the mass flux of dissolved CO 2 passing through the elements along the A -A' profile at specified times. Figures 15 and 16 show comparisons of dissolved CO 2 mass fluxes for the test cases. We focused on the dissolved CO 2 mass fluxes after 30 days to exclude the initial diffusion period (Fig. 15). According to previous investigations this diffusion period might influence the convection fluxes calculation along the A -A' profile. Figure 15a depicts the temporal variations of average CO 2 mass flux cross the A -A' profile, and Fig. 15b accumulates the dissolved CO 2 mass passing the profile. The average CO 2 mass fluxes decreased gradually and ultimately tended to reach stable values later. The results in Fig. 15a clearly show that an increase in lnk 2 v can lead to an increase in average mass fluxes. Different corre-lation ratios (i.e., Cases 5, 6, and 7) yielded relatively small average and total mass fluxes according to a comparison of the calculated mass flux and total mass in Cases 5, 6, and 7 with Case 2. Detailed quantifications of the differences are presented in Fig. 16.
In Fig. 16 the bar charts show comparisons of the accumulated total mass of dissolved CO 2 that passes through the A -A' profile for 30 and 200 days. When Cases 3 and 4 were compared with Case 2 the results clearly show that permeability variations can enhance CO 2 convection. Our simulation results indicate that an increase in lnk 2 v from 1.5 (Case 3) to 2.0 (Case 4) can slightly enhance CO 2 convection from 2.5% (Case 3) to 5.6% (Case 4) in comparison with that in Case 2. The CO 2 cumulative mass increased from 3.5 -8.1% after 200 days. Increasing the lateral small-scale lnk correlation slightly reduced the vertical CO 2 convection activity according to a comparison of the results (i.e., Cases 5, 6, and 7) with those in Case 2 (negative percentages for Cases 5, 6, and 7).
Notable CO 2 convection was obtained in Case 8 because the mean lnk was one order of magnitude greater than that in Case 2. The mean lnk value for Case 9 was one order of magnitude lower than that for Case 2. After 30 days the cumulative CO 2 mass for Case 9 was 12.7% smaller than that for Case 2. This difference decreased to 50.3% after 200 days. Cases 8 and 9 used the same statistical structure to generate RFs. The behaviors are different in both fingering patterns (Fig. 14) and the convection flux magnitudes (Fig. 16). For a specified location the dissolved CO 2 mass flux was influenced by local small-scale variability in permeability. A positive sign indicates that the masses of the cases were greater than that of Case 2. However, a negative sign represents that the masses of the cases were less than that of Case 2. This permeability contrast triggered fingerings with short durations. In addition to the variations in small-scale permeability the mean permeability of a reservoir formation is a crucial parameter for controlling long-term convection.

COnCluSIOnS
We propose one-way downscaling processes for conducting CO 2 migration numerical simulations at different scales. The hierarchy tree in the downscaling processes enabled us to extract solutions from the large-scale main model and a series of sub-models. We developed a data integration program for formulating the input data from the TOUGH-REACT model for each specified sub-model at specified times. Without global grid refinements for small-scale problems the simulations can directly focus on the interested sub-models. The downscaling processes are computationally efficient and not limited to the TOUGHREACT model used in this study for flow and transport simulations.
Our numerical examples showed that fingerings always occur below the bottom edge of the transition zone. This is crucial information for facilitating quantification of the additional area in dissolved CO 2 convection at various degrees of permeability variations and times. The initiation and development of fingerings depend highly on the instability triggered by local permeability variations. An increase in lnk variance leads to an increase in CO 2 convection compared with the homogeneous case. In each numerical case CO 2 convection was relatively active in the early time periods and ultimately reached a stable value in late time periods.
Numerical cases for different anisotropic ratios revealed that an increase in lateral correlation for permeability variations might reduce the dissolved CO 2 convection. However, a decrease in CO 2 mass flux with an increase in anisotropic ratios was not significant. The local variations in permeability near the lower XCO2a transition zone boundary control fingering occurrences and development. At the same lnk variance the RF generator produced lnk fields with highly lateral correlations and this might yield slightly smaller permeability variations near the bottom edge of the transition zone. Therefore, anisotropic cases produced relatively small CO 2 convections in this study. The mean permeability plays a crucial role in transporting dissolved CO 2 . Identical permeability variations in a small mean permeability reservoir can lead to a delay in fingering occurrences and a small CO 2 mass flux. In the numerical experiments conducted in this study, the permeability variations controlled convection in the early time periods. In the long convection term the effects of average permeability in a reservoir became the significant properties.
We focused on the dissolution trappings of CO 2 migrations, particularly regarding the influence of permeability variations in reservoir materials. Moreover, other factors involved in mechanical and chemical interactions can produce instabilities in dissolution trappings. In addition to permeability variations, the mechanical and chemical interactions might improve or deactivate the development of fingerings. The developed downscaling framework enables conducting further investigations into the effects of other factors on CO 2 convection.