Effect of permeability heterogeneity on the dissolution process during carbon dioxide sequestration in saline aquifers: two-and three-dimensional structures

The convection–diffusion process of carbon dioxide (CO2) dissolution in a saline reservoir is investigated to shed light on the effects of the permeability heterogeneity. Using sequential Gaussian simulation method, random permeability fields in two and three-dimension (2D and 3D) structures are generated. Quantitative (average amount of the dissolved CO2 and dissolution flux) and qualitative (pattern of the dissolved CO2 and velocity streamlines) measurements are used to investigate the results. A 3D structure shows a slightly higher dissolution flux than a 2D structure in the homogeneous condition. Results in the random permeability fields in 2D indicates an increase in the standard deviation of the permeability nodes enhances the dissolution efficiency, fluctuations in CO2 dissolution flux, separation between the different realizations from the same input parameters, and tendency toward more jagged convective fingers’ shape. Furthermore, the distance between the permeability nodes increases the convective fingers’ dissolution efficiency and jagged structure. The degree of freedom in 3D structures results in a higher chance of escaping from the low permeability zones and reduces the interactions between convective fingers in 3D systems. With the same variance and correlation length between permeability nodes, connectivity between high permeable zones in 3D cases are less than that of 2D cases; therefore, 2D realizations overestimate the dissolution flux of real heterogeneous 3D structures, which should be considered carefully. CO2 sequestration in two and three dimensional heterogeneous saline aquifers are investigated. 3D structures in homogeneous conditions show higher dissolution than 2D structures. 2D realizations overestimates the dissolution flux over real heterogeneous 3D reservoirs. CO2 sequestration in two and three dimensional heterogeneous saline aquifers are investigated. 3D structures in homogeneous conditions show higher dissolution than 2D structures. 2D realizations overestimates the dissolution flux over real heterogeneous 3D reservoirs.

into the Earth's atmosphere (Dai et al. 2016(Dai et al. , 2020Singh et al. 2019Singh et al. , 2020Sturmer et al. 2020;Wang et al. 2020). CO 2 dissolution in saline aquifers promotes natural immobilization of CO 2 due to the specific properties of saline aquifers such as chemistry Fu et al. 2020;Soltanian et al. 2019), porosity, temperature and pressure (Amarasinghe et al. 2020), and huge CO 2 storage potential (Ershadnia et al. , 2021Li and Jiang 2020;Mahmoodpour and Rostami 2017;Mahmoodpour et al. 2018;Tang et al. 2019). The CO 2 injected into saline aquifers is typically 40-60% lighter than reservoir brine (Lemmon et al. 2019). Due to buoyancy forces, the CO 2 plume moves upward and spreads under the cap-rock. The leakage of low density and low viscosity CO 2 through the possible permeable zone, is a concern for the project's success. The dissolution into underlying brine and subsequent mineral immobilization helps the security of the storage site (Islam et al. 2013;Emami et al. 2015a). Mineral trapping happens in a very long time-scale in comparison with dissolution, especially in sandstones (Moghaddam et al. 2012;Celia et al. 2015;Omrani et al. 2021). Thus, dissolution trapping has a critical role in mitigating anthropogenic CO 2 emission (Jia et al. 2020). CO 2 dissolution in brine not only enhances the aquifer capacity for CO 2 storage but also reduces the geomechanical risks with respect to the storage integrity, thus helping the security of the project (Kong and Saar 2013).
In the early period, injected gas dissolves into the pore fluid through the diffusion mechanism and becomes denser than the original brine. The presence of denser fluid over the reservoir fluid initiates gravitational convective instabilities (Riaz and Cinar 2014;Hassanzadeh et al. 2005;Farajzadeh et al. 2007a). In proper conditions, instabilities evolve and generate large-scale convective flow. Convective mixing appears in the downward-traveling fingers . Slumping of the dense gravity current controls the dissolution rate because it sweeps fresh water upward causing increase in the concentration difference at CO 2 -brine interface (Macminn et al. 2012). The convective mixing process is highly unstable due to the continuous formation and coalescence of fingers (Xie et al. 2011), and traditional reservoir simulation tools cannot handle the small-scale trapping dynamics at the geologic basin scale (Hidalgo et al. 2013).
It is fair to say that most of the studies that simulate CO 2 dissolution in saline aquifers have used 2D homogeneous models. Previous studies examined the numerical modeling of convective mass transfer (Farajzadeh et al. 2009;Mahmoodpour et al. 2019a), convective instability (Slim 2014;Emami et al. 2015b), impurities (Li and Jiang 2014;Li et al. 2015;Mahmoodpour et al. 2020;Soltanian et al. 2018) as well as experiments on density-driven convection (Kneafsey and Pruess 2010;Mojtaba et al. 2014;Khosrokhavar et al. 2014;Mahmoodpour et al. 2019a), and Hele-Shaw flows de Paoli et al. 2020). Despite the more realistic nature of heterogeneous media, these have received less attention. Time-dependent velocity oscillation due to permeability variation leads to different dissolution patterns in heterogeneous systems (Islam et al. 2014). Reservoir heterogeneity may arise from two aspects: structural heterogeneity due to faults and folds or salt diapers and heterogeneous stratigraphic facies (Lengler et al. 2010;Liu et al. 2020). Furthermore, reducing the dimension of study from 3D to 2D could help in examining the field scale problem. Nevertheless, the effect of this dimension reduction is unclear because of limited number of 3D studies focused on structural heterogeneity.
Regarding 2D simulations, Green and Ennis-King (2010) used small-scale barriers to create heterogeneous porous media. They found that it is possible to approximate the heterogeneous reservoir with an equivalent anisotropic homogeneous reservoir. Lengler et al. (2010) developed a heterogeneous aquifer model for Ketzin with the variance of 0.55-1 for the permeability range of 0.02-5000 md to examine the injection stage. They concluded that a heterogeneous porous medium gives higher storage capacity than the homogeneous one. Ranganathan et al. (2012) produced permeability fields with a standard deviation of 0.1-2 and correlation length of 0.025-0.1 m. They observed three dissolution regimes as fingering, channeling, and dispersive based on the standard deviation and correlation length values. Kong and Saar (2013) observed an increasing trend of dissolution flux with permeability heterogeneity. Chen et al. (2013) concluded that increasing heterogeneity retards the onset of convection because low-velocity regions resulting from low-permeability zones decrease the system's instability. Furthermore, they reported that the variation of permeability zones has a controlling role in the gravitational process and correlation length has a negligible effect. However, the dissolution rate is higher in heterogeneous media, and correlation length significantly affects dissolution rate at later times (Chen et al. 2013). Elenius and Gasda (2013) extended the work of Green and Ennis-King (2010) to large-scale barriers resulting from mudstone inclusions or shale deposits. They showed that mass flux decreases with barrier length and increases with distances between them. Furthermore, it is possible to estimate the dissolved CO 2 flux in these configurations too (Elenius and Gasda 2013). Islam et al. (2014) simulated the nonisothermal dissolution of CO 2 in in a saline aquifer and pointed out that heterogeneous porous media has higher CO 2 dissolution flux than homogeneous one. Agartan et al. (2015) examined the effect of heterogeneity experimentally using similar fluids and concluded that heterogeneous formations may not always enhance the convective mixing. Singh et al. (2021) found that heterogeneity in permeability may enhance the overall CO 2 dissolution in brine. Regarding 3D systems, researchers mainly focused on the structure-based heterogeneity. They showed that the low permeability regions might decrease the vertical spreading, but they may enhance transverse CO 2 plume spread Ershadnia et al. 2020). Furthermore, there are limited 3D studies in homogeneous systems (Pau et al. 2010;Li and Jiang 2014). Therefore, there has been ongoing ambiguity regarding CO 2 dissolution patterns in complicated 3D heterogeneous structures.
To address these open issues, we performed the numerical study. We examined whether it is possible to estimate the behavior of the heterogeneous systems in 3D and changes in the behavior of dissolution flux resulting from dimension reduction from 3D to 2D.
The following structure is adopted for this paper: the methodology describing the governing equations, selected domain, boundary conditions, and methodology for generating heterogeneous media. In the results and discussion section, first, the ideal structure of permeability results are discussed. An ideal system's results help us explain the different behavior in random fields. We then examine the random permeability fields, focusing on the effects of permeability distribution variance and distance between permeability nodes on the dynamics and mixing efficiency of the process in 2D systems. We then extend our study to a more complicated 3D system.

Methodology
A rectangular cubic space far away from an injection well, where the horizontal gas-brine contact is selected to be the domain of interest to capture the actual behavior of the system. We assumed a height of 10 m (Bachu et al. 2005) and a length of 20 m for our system. It should be noted that this value is much higher than the diffusive boundary layer thickness for these conditions (Kong and Saar 2013). There is no CO 2 present initially in the system and the top boundary is assigned with a constant concentration of 1 where the maximum concentration is measured relative to the initial solubility equilibrium. A sharp interface is considered for the CO 2 -brine contact. This is a valid assumption for reservoirs deeper than 1 km (Ennis-King et al. 2005;Islam et al. 2014). It means we assume that there is no transition zone of saturation of brine and CO 2 phases in the top section of the system (no capillary effect), and we consider only flow and transport in the aqueous phase. We considered that flow is possible only through the top boundary. The no-flow boundary condition ensures that the simulation domain shows a laterally infinite system. Similar boundary conditions were implemented in previous studies (Ranganathan et al. 2012;Islam et al. 2014;Rezk and Foroozesh 2019). Furthermore, Ranganathan et al. (2012) examined different side and bottom boundary conditions effects on the process. They concluded that dissolution flux with different boundary conditions has the same trend (we have shown the same result in Sect. 3.3 and Fig. 15). However, at late times, the periodic boundary condition for lateral sides may show slightly higher dissolution flux. Chemical reaction simulations are mainly relevant for long-term geological CO 2 sequestration investigations (Dai et al. 2020); therefore, in the considered timescale of this study, we can neglect the effect of chemical reactions. The following equations are coupled and solved using finite element method in COMSOL Multiphysics: Mass transfer equation: (1) where ∅ , C (mol/m 3 ), D , and u show porosity, dissolved gas concentration, molecular diffusion coefficient, and the Darcy velocity vector, respectively. Fluid flow: where k , , P , , g , and H are permeability, viscosity, pressure, density, gravity acceleration constant, and height respectively. The density of brine plus dissolved gas is calculated with Eq. 3 (Garcia 2001): subscripts aq , w , and i indicates solution, water, and solute. Also, M , x , NCG , and V represent molar weight, solubility values, number of dissolved components, and partial molar volumes. Table 1 presents the base case data which are used during the simulation. It should be noted that small flow velocities reduce hydrodynamic dispersion (D) to the molecular diffusion (Ranganathan et al. 2012). Parameters are selected based on reported data from previous studies to have a representative model from the actual conditions. For example, considering permeability values, Law and Bachu implemented permeability values in the range of 100-400 md (Law and Bachu 1996), Lengler et al. (2010) used an average value of 90 md and  built permeability fields in the range of 500-5000 md. Maximum CO 2 solubility in brine is calculated from a thermodynamic model based on pressure and temperature values (Hassanzadeh et al. 2008). Dissolved concentrations is expressed in the range of 0-1 relative to this maximum solubility. In real porous media, heterogeneity in permeability, porosity, or fluid properties initiates instabilities, but in ideal homogeneous cases, we must impose instabilities on the system. For this purpose, we follow one of the most common schemes, and thus, sinusoidal perturbation is implemented on the top boundary. The long-term behavior of a system is independent of this perturbation (Farajzadeh et al. 2007b).
The Rayleigh number is an important factor in analyzing convection-diffusion problems (Lapwood 1948): In which k , Δ , g , H , , ∅ , and D are permeability, density difference, gravity acceleration, height, viscosity, porosity, and the diffusion coefficient, respectively. Sequential Gaussian Simulation (SGS) is a common algorithm to create a heterogeneous porous medium (Delbari et al. 2009;Safikhani et al. 2017;Jia et al. 2018). We employed the SGS method through SGeMS software to create a random permeability field (Bianchi and Zheng 2009). SGeMS generates normal random distribution in the range of 0-1. Because permeability distribution is commonly lognormal, Eq. 6 is used to create log-normal permeability distribution from the normal distribution resulting from SGeMS (Nield et al. 2010;Ranganathan et al. 2012): Here, k log −normal , k average , and k SGeMS shows permeability field in log-normal distribution, average permeability, the standard deviation, and correlation length obtained from SGeMS. According to a specified correlation length, these random permeability values are assigned to the simulation domain. The permeability of space between nodes is interpolated to result in a continuous field. To examine the effect of the variance in the permeability field, permeability values with the standard deviation of 0.25, 0.5, and 1 are distributed in 2D space with 1 m distance in x and y directions. Furthermore, to examine the effects of the distance between permeability nodes, three sets of simulations with a correlation length of 0.5, 1, and 2 m and = 1 are compared. While permeability fields are generated randomly, therefore; for each set of simulations 20 realizations are considered.
(4) It should be noted that mean permeability is the same as reported in Table 1. To compare results between 2D and 3D cases, 10 realizations in a 3D case with the standard deviation of 1 and a distance of 2 m are generated. Furthermore, 15 ideal cases are considered to examine the effect of high-permeability zones locations and their connectivity on the dissolution process (Fig. 2). The results of these ideal cases will be used to explain the behavior of random fields. In Fig. 2, light regions show permeability of 200 md and dark ones show permeability of 400 md. It should be noted that previously this methodology is used to simulate experimental tests results with the authors (Mahmoodpour et al. 2020).

Results and discussions
Two Measurements involving the average concentration of dissolved CO 2 and dissolution rate, with visual dissolution patterns, are used to examine the process. The initial concentration of CO 2 in water is assumed to be zero. We used a thermodynamic model based on temperature and pressure to calculate the CO 2 solubility in water. We normalize the dissolved CO 2 concentration based on maximum solubility for the given P − T conditions to make it vary from zero to one. Areal or volumetric average gives dissolved CO 2 mass in the normalized form. Then dividing this value by the surface area or volume gives the average concentration of dissolved CO 2 . We calculate the dissolution flux by using the difference of dissolved CO 2 between two consecutive time steps.
Previously, we used Rayleigh number to estimate the shutdown regime based on the experimental results (Mahmoodpour et al. 2019b). The Rayleigh number of our system in this work is around 2000 that gives us the time of 5.2 years for the start of the shut-down regime (Mahmoodpour et al. 2019b). Therefore, simulations are followed for up to 7 years when the convective front hits the system's bottom. The proposed methodology is validated based on the experimental data (Mahmoodpour et al. 2019b(Mahmoodpour et al. , 2020. Base case simulation results are provided in Fig. 1, for which inputs are taken from Table 1. Figure 1a shows the dissolution pattern and velocity streamlines at 3 years in which convective fingers are equally spaced. Because the lateral boundaries are considered as the no-flow boundary conditions, when convective fingers move laterally and reach these  Table 1: a dissolved CO 2 pattern at 3 years; b dissolved CO 2 pattern at 7 years; c mass averaged dissolved CO 2 ; d dissolution rate boundaries, their movements in the lateral direction are restricted and force them toward the bottom. However, for an infinite system, this may not be true, and therefore, our model slightly increases the dissolved CO 2 mass (In Sect. 3.3, we have shown that this increase is marginal). Therefore, these fingers grow more in comparison with other fingers As time progresses, tiny fingers coalesce to each other and create larger and stronger (with higher velocity toward the bottom of the system) ones. This behavior is depicted in Fig. 1b in which the number of fingers decreased at 7 years. Figure 1c shows the average value of the dissolved CO 2 . At early times, the averaged dissolved CO 2 grows at a higher rate than the later times. This behavior is more clear in Fig. 1d, where depicts the dissolution flux. As time goes on, the concentration gradient driving force to the diffusion mechanism decreases, and dissolution flux reduces. However, during specific periods, activation of the convective mechanism tries to prevent dissolution flux decrease, as depicted in Fig. 1d. To detect these regimes in the dissolution flux curves, a sharp declination (diffusion Fig. 2 Schematic of the high permeability zones distribution in ideal cases. Light regions shows permeability of 200 md and dark region represents permeability of 400 md dominated regime) is followed by an increase of dissolution flux (flux growth regime), a plateau (quasisteady-state regime), and finally, a gradual decrease of dissolution flux (shut-down regime) (Mahmoodpour et al. 2019a). In the following sub-sections, first, the system's behavior in ideal heterogeneous porous media is examined. Findings from this part are used to analyze the system's behavior in more complex heterogeneous porous media in 2D and 3D cases.

2D ideal cases
Simple cases are used to analyze the impact of a high-permeability zone's structure on the dissolution efficiency. In addition, ideal 2D cases may be used to explain the differences between more realistic random fields. Our aim is to gain insight into a high-permeability zone's location and connectivity impact on CO 2 distribution. Figure 2 represents 15 cases that are used in this regard. Darker blocks indicate a permeability of 400 md. The first column shows schematics used to examine the layering impact. The no-flow boundary condition near lateral walls accelerates finger movement in this region and slumping of dissolved CO 2 restricts finger migration upward. Therefore, only downward fingering is possible. To study these arrangements, the second column has the same quantity of high permeability blocks in a vertical direction. The third column is employed to examine the connectivity effect.
The average amount of the dissolved CO 2 and dissolution flux for these 15 cases are represented in Fig. 3a, b, respectively. Convective fingers tend to slump through high permeability zones; therefore, more interaction in such regions. A highly permeable layer in the first row of case 1 results in a faster merging of small fingers and the highest dissolution rate at early times. Lowering the highly permeable layer location hampers the process. Based on this rule, case 5 shows the lowest efficiency despite its good connectivity. This behavior has been confirmed visually with experimental studies, in which the presence of high permeability zones at the top of the system enhanced the dissolution process (Agartan et al. 2015).
The second column in Fig. 2 reveals the locations of the highly permeable columns through the system. In the proximity of the sidewalls, finger movement is constricted to the downward direction; therefore, we expect the highest dissolution in this situation, which is confirmed by the results. By receding of high permeability zones from the sidewalls, the dissolution capacity reduces until the distances of these two highly conductive pathways reach a critical point. It is reasonable to imagine that dissolution is enhanced with a large number of surviving fingers. However, surprisingly, cases 9 and 10 show different behaviors. These cases reveal the importance of the connectivity of high permeability zones to allow the merging of tiny fingers. Based on this explanation, the structure of case 10, in contrast to the most significant distance from the sidewalls, has the most favorable efficiency through cases 6-10.
The third column of Fig. 2 highlights the impact of connectivity. In the schematic of case 11, despite locating five blocks in the first row (which causes the same onset time as case 1 for convective flow), the amount of dissolved gas in this situation is smaller than the average of cases 1 and 2. In another exciting configuration, case 14 uses the advantage of being near the first row and sidewalls. Considering all cases, the flux growth regime (see Mahmoodpour et al. 2019b) is more evident in case 1 and case 14. However, the dissolution flux of case 1 decreases rapidly because the permeability of the lower parts of the system is small. Case 14 is a combination of the conditions which may result in an efficient dissolution process. In this case, high permeability zones are located at the top surface or close to the side boundaries and there is a good connectivity between these zones. But in case 1, these high permeability zones are only located at the top section, which is the controlling item for merging the convective fingers at early times. For these reasons, at the later times, the dissolution efficiency of the case 14 is higher than in case 1. Pattern of the dissolved gas for the cases with the maximum and minimum dissolution efficiencies are shown in Fig. 3c-h. The results reveal that the spatial distribution of high permeability zones and their correlation lengths are important factors to analyze CO 2 dissolution behavior in heterogeneous structures. We then examined more realistic random field systems based on these insights from ideal structures.

2D Random field
Random permeability fields are used to verify heterogeneity effects in more realistic cases. Three levels of standard deviation (0.25, 0.5, and 1) and distance (0.5, 1, and 2 m) are considered to examine the effects of standard deviation and distance between permeability nodes. The permeability node represents permeability assigned at that node and averaged over a zone specified by correlation length. Because the distribution pattern from a realization may be very different from another one with the same inputs, it is essential to repeat the simulations in each case for many realizations to be sure about the conclusion. As this is not practically feasible, we have used 20 realizations for each case due to the high calculation time. Figure 4a, b demonstrate the average amount of the dissolved CO 2 (normalized based on the CO 2 solubility in water) and dissolution flux curves for realizations with the standard deviation of 0.25 and distance 1 m. Results reveal that all realizations follow up a similar trend of fast decline of dissolution flux at early times and then activation of the convective flow, which keeps the dissolution flux constant. The average amount of the dissolved CO 2 after 7 years is lower than 0.25 of CO 2 solubility at this pressure and temperature for all realizations. The separation between the dissolution flux curves is slight, and they are distributed in a slim band. Realizations with minimum and maximum amount of the dissolved CO 2 are examined in Fig. 4c-h. Comparison between the permeability field of the minimum dissolved CO 2 (Fig. 4c) and the permeability field of the maximum dissolved CO 2 (Fig. 4d) shows that the presence of the higher permeability nodes in the top section of the reservoir, which is close to the gas-liquid interface is responsible for this dissolution flux separation. This behavior is observed in ideal cases too. A comparison of the dissolved CO 2 patterns in 7 years shows that the convective fingers move faster in the maximum case (Fig. 4f) compared to the same time of the minimum case (Fig. 4e). Furthermore, the interaction between convective fingers is higher in the maximum case, creating giant convective fingers. This phenomenon is more at late times in which, after 7 years of the process, there are a few large convective fingers in the maximum case (Fig. 4h) compared to the minimum case (Fig. 4g). Because of the slight standard deviation between the permeability nodes (permeability range is 50-450 md in these Figures), the convective fingers do not move toward the high permeability zones and move almost vertically toward the bottom of the system. This results in nearly symmetric velocity streamlines.
The average amount of the dissolved CO 2 and dissolution flux curves for realizations with the standard deviation of 0.5 and correlation length of 1 m is shown in Fig. 5a, b. Curves show higher separation from each other compared to the realizations with a standard deviation of 0.25. Furthermore, some cases reach to 0.25 of the mass averaged dissolved CO 2 . In addition to this observation, the total dissolved mass in the minimum case is higher for the standard deviation 0.5 compared to the standard deviation of 0.25. However, the difference between these minimum cases is small (0.2099 for the standard deviation of 0.25 and 0.2139 for the standard deviation of 0.5). These observations show that increasing variance (intensifying the difference between permeability nodes) increases the dissolution rate for the investigated cases. Examining the dissolved CO 2 pattern in Fig. 5e-h reveals that the convective fingers at early times in the maximum realization of the standard deviation of 0.5 (Fig. 5f) is smaller than the related pattern the standard deviation of 0.25 (Fig. 4f). This emphasizes the effect of high permeability nodes on creating the strong convective fingers at early time, which finally enhances the final amount of the dissolved CO 2 . The creation of a strong finger at early time in Fig. 5h results in other convective fingers not developing as much as Fig. 4h. This strong finger movement in front reduces the concentration gradient in the vertical direction for other fingers that drives the convective finger development. Furthermore, due to the higher contrast between permeability nodes compared to the realizations with a standard deviation of 0.25, there is a slight tendency of the most enormous finger's movement toward the high permeability zones. Figure 6a, b show the average amount of the dissolved CO 2 and dissolution flux curves for realizations with the standard deviation of 1 and a distance of 1 m. Differences between the realizations with the Vol.: (0123456789) same distance and different standard deviations discussed for the standard deviation of 0.25 and 0.5 are intensified by increasing the standard deviation to 1. The maximum value for the average amount of the dissolved CO 2 approaches to near 0.3 and dissolution fluxes separate from each other more than before. Furthermore, the final amount of the dissolved CO 2 for the minimum realization with the standard deviation of 1 is 0.2201, which confirms the previous trend of the increasing dissolution flux with increasing the variance. Dissolved CO 2 patterns in Fig. 6f, h show more jagged shapes compared to the related patterns in Figs. 4 and 5 resulting from a high difference between permeability nodes. As a result of these high differences between permeability nodes, the tendency toward the high permeability zones increases and creates these jagged shapes. This tendency is more in the bottom right corner of Fig. 6h. High permeability nodes accumulate and the convective finger deviates from its' vertical movement toward the system's bottom. This deviation in fingers' movement makes the velocity streamlines more asymmetric compared to cases with lower variance values. In this study, the porosity values are kept constant in all parts of the system. Consider a scenario where we do not consider constant porosity values and use a experimentally measured k − dependency. In that case, the preferential flow toward high permeability zones will be decreased. Because considering Rayleigh number (Eq. 5) compares the effectiveness of the convective flow to the diffusion flow and controls fluid movement toward the high permeability zones. In this equation, permeability and porosity have diverse impacts and follow the same increasing or decreasing trend. Therefore, considering a relationship between porosity and permeability will reduce the fluctuation of the Rayleigh number in different cases. However, conclusions would be the same as this study; but the degree of differences between cases would be decreased.
To examine the impact of the correlation length, only the of 1 is considered. Figure 7a, b demonstrate that average amount of the dissolved CO 2 and dissolution flux curves for realizations with the standard deviation of 1 and the correlation length of 0.5 m. Compared with correlation length 1 m, there is a higher separation between the dissolution flux curves. This behavior results from a large number of the permeability nodes in correlation length 0.5 m, where high permeability nodes increase the dissolution flux and low permeability nodes decrease this value. However, the final amount of the dissolved CO 2 is lower than realizations with 1 m correlation length. Because the high permeability nodes in 0.5 m correlation length case contribute to the small space of the system, they are surrounded by low permeability zones. This effect is more evident in dissolution patterns depicted in Fig. 7e-h. The comparison with the related patterns in Fig. 6 reveals that there are weak interactions between convective fingers at early times resulting in a high number of convective fingers.
Furthermore, convective fingers have smooth structures and move vertically toward the bottom of the system with a lower deviation toward the high permeability zones. These patterns and dissolution amount show the connectivity between high permeability zones, which are explained in ideal cases. With realizations in the distance 0.5 m, high permeability zones are not well connected or concentrated compared to those of length 1 m. This poor connectivity eventuates in the lower amount of the dissolved CO 2 and well smoothed convective fingers pattern. To go deeper into the effect of the permeability nodes distance, 20 realizations with a length of 2 m are created and the simulation process is repeated with them. Figure 8 represents the quantitative and qualitative measures of this process. Figure 8a shows that realizations with a distance of 2 m have a greater averaged dissolved mass compared to 1 m correlation length (Fig. 6a). A possible reason for this difference is explained in the last part about the connectivity of the permeability nodes and captured during ideal case analysis. Furthermore, due to the localised heterogeneity, the fluctuation between dissolution fluxes is more in comparison with the distance of 1 m. This behavior is more in Fig. 8c-h, which shows the dissolution patterns where convective fingers have more jagged shapes and wider horizontal scatter compared to the distance of 1 m (Fig. 6c-h). Strong interactions between the convective fingers result in some strong convective fingers at early times. The tendency toward the high permeability zones is stronger, resulting in a wider plume pattern in the horizontal direction and shaping more asymmetric velocity streamlines compared to the distance of 1 m.

3D structures
To examine the difference between 2 and 3D structures, first, simulation is done with base case input in a homogeneous condition (Table 1). Quantitative measures of the dissolution process are depicted in Fig. 9a, b. Even the overall behavior and trends of the curves are similar at early times, the dissolution process in the 3D homogeneous case is slightly higher than the 2D homogeneous case. Similar behavior is observed with Pau et al. (2010), where they examined the CO 2 dissolution in a homogeneous porous media. They concluded that more degrees of freedom in 3D structures result in a more complex fingering phenomenon. Especially when convective fingers reach the top of the system and spread laterally. In the 3D design, they move in all directions compared to the 2D system, where they move along a horizontal line. For 3D cases, some narrow ridges with higher downward velocity develop at the fluid interface. This is due to the higher dissolution flux than the 2D cases. Figure 9 shows the dissolution patterns in a 3D homogeneous system. Figure 9c, d reveal that different cross-sections from the 3D results have different dissolution patterns at early-and late times; therefore, the 2D and 3D dissolution pattern do not match. Isosurface plots of the dissolved CO 2 concentration (Fig. 9e, f) show that convective fingers have almost similar shapes and are distributed symmetrically in the system except for the fingers that are close to the boundaries. Convective fingers horizontal movement restriction resulting from the no-flow boundary condition is more evident in velocity streamlines in Fig. 9g. The finger created in the corner of the system experiences constraints from two sides walls; therefore, its' movement toward the bottom of the system intensifies. However, inner convective fingers develop similar velocity fields around themselves. At late times where all convective fingers are well developed, the similarity between streamlines is increased (Fig. 9h). Velocity values toward the system's bottom are higher than the upward movement because the volume of down-ward moving convective fingers is smaller than the upcoming fluids. Furthermore, velocity values at the tip of the convective fingers are higher than the finger's root resulting from the conical shape of the convective fingers.
Due to the high computational time of the simulations in 3D heterogeneous porous media, ten realizations with the standard deviation of 1 and distance of 2 m are considered. The average amount of the dissolved CO 2 and dissolution flux curves are shown in Fig. 10. 3D realizations have lower fluctuations in dissolution flux than 2D cases. This is because in 3D structures, high permeability zones are more surrounded by low permeability zones, and low permeability zones are more surrounded by high permeability zones than similar realizations in 2D cases. In other words, with the same variance and correlation length of the permeability nodes as the inputs for SGS, the connectivity between the high permeability zones in 3D cases are lower than in 2D cases. Especially at early-and late-time dissolution rates show very close behavior. Because, in these periods, the flux is mainly controlled with diffusion coefficient (as a function of T and P). Based on the finding results from the ideal cases, lower connectivity in 3D eventuates in the lower amount of the dissolved CO 2 in 3D cases. From the qualitative analysis of the realizations, the permeability field of the realization with the lowest-and the highest amount of the dissolved CO 2 are depicted in Figs. 11a, b for top views and Fig. 11c, d for the bottom views. As it is expected from the analysis in previous parts, in the realization with the highest amount of the dissolved CO 2 , the high permeability zones are concentrated in the top section of the system (close to the gas-liquid interface). This configuration helps to interact with the initial weak convective fingers to create stronger ones at early times at high permeability zones. Cross-sections of the dissolved CO 2 patterns in Fig. 12 show this phenomenon. Figure 12b represents that the positions of the strong convective fingers are well aligned with the positions of the high permeability zones, especially for the zones close to the no-flow boundaries. These zones use the advantages of the high permeability and lateral movement restriction; therefore, the strongest convective finger develops.
Isosurface plots for dissolved CO 2 patterns in Fig. 13 show better convective fingers' development in high permeable zones. As it is expected, the convective fingers are not as equally distributed as the homogeneous cases and they show some jagged structures similar to the related 2D cases. However, higher degrees of freedom in 3D structure smoothens the convective finger's shape. The velocity streamline in Fig. 14 shows that convective fingers in high permeability zones and close to the no-flow boundaries have higher velocity than other fingers.
In this section, three different boundary conditions are compared (a) No-flow condition for bottom and side boundaries, (b) Periodic side boundary conditions and impermeable bottom boundary, and (c) Zero concentration at the bottom boundary with no-flow at the side boundaries. Figures 15a, b compare three different boundary conditions for variance 1 and correlation length 1 and 2 m, respectively. It is evident that the average amount of dissolved CO 2 is approximately the same for all three boundary conditions. However, at later times, the no-flow boundary condition shows slightly higher dissolved CO 2 probably due to pressure buildup.

Conclusions
Analysis of the CO 2 dissolution in ideal heterogeneous cases reveals that position and connectivity of the high permeability zones are two essential factors that control the amount and pattern of the dissolved CO 2 in porous media. Increasing the variance of the permeability nodes in 2D structures with the same distance between the permeability nodes increases the dissolution efficiency, fluctuations in dissolution flux, and separation between the different realizations from the same input for the permeability field. Furthermore, with increasing the difference between permeability nodes (variance), the tendency toward the high permeability zones increases, and it changes the shape of the convective fingers from the smoothed to the jagged ones. It deviates their movement direction from vertical toward high permeability zones. Distance between the permeability nodes increases the dissolution efficiency and irregular structure of the convective fingers too. But, a general role for the relationship between the dissolution flux fluctuation amount and distance between nodes is not observed.
Considering the homogeneous permeability fields, the dissolution process in 3D homogeneous cases is slightly higher than in 2D homogeneous cases. More degrees of freedom in 3D structures result in a more complex fingering phenomenon and formation of some narrow ridges with higher downward velocity than the 2D system. Therefore, 2D simulations may underestimate the dissolution flux of 3D homogeneous structures. In heterogeneous realizations, 3D realization shows lower fluctuations in dissolution flux and makes their behavior predictability easier than 2D cases. If we do not consider constant porosity values and use a measured relationship between the permeability and porosity from core analysis to assign porosity based on the permeability field, the preferential flow toward high permeability zones will be decreased and fluctuations may become smaller. Therefore, using a homogeneous porosity field enables us to use the limited realizations in 3D to predict the dissolution flux.
Furthermore, lower connectivity between high permeability zones in 3D structures results in lower dissolution flux than 2D structures. Therefore, this overestimation in dissolution flux from 2D simulations in heterogeneous porous media should be considered carefully. Owing to the CO 2 dissolution, the overall pressure of the system and consequently the possibility of leakage decreases. Therefore, using the results of 2D realizations may eventuate in some over-pressurized periods during CO 2 sequestration projects in heterogeneous porous media.