The governing equations for convective mixing include the continuity equation and the mass balance equation. The continuity equation in the porous medium is given by ∂ ϕ m ρ ∂ t = − ∇ ∙ ρ u (1)

Darcy’s velocity u is calculated by Darcy’s law u = − K m μ ∇ p − ρ g ∇ D (2)

Where t is the time, ϕ m is the matrix porosity, ρ is the mixture density, K m is the matrix permeability, μ is the dynamic viscosity of the mixture, p is the pressure and g is the gravity term. The mass balance equation for the transport of dissolved CO₂ and H₂S is ∂ ϕ m c i ∂ t = ∇ ∙ ϕ m D e i ∇ c i − c i u , i = C O 2   o r   H 2 S (3)

Where c i is the component concentration and D e i is the effective diffusion coefficient of the component. The Millington and Quirk model (Millington and Quirk, 1961) is used to calculate the effective diffusion coefficient given by D e i = ϕ m 4 / 3 D i (4)

Fluid density is a function of temperature, pressure, and component concentration. The density equation expressed by Diersch and Kolditz (2002) is used for simplicity. ρ = ρ 0 1 + ∑ i = 1 n β i c i − c 0 , i (5) β i = 1 ρ 0 ∂ ρ ∂ c i T (6)

Where ρ 0 is the density of CO₂-free water, β i is the densification coefficient.

To account for the presence of fractures within the porous medium, the “Fracture Flow” interface available in COMSOL Multiphysics 5.6 was utilized. This interface utilizes tangential derivatives to characterize the flow along the interior boundaries. The flow in the fracture is governed by the continuity equation and the mass balance equation, represented as (COMSOL, 2020). d f ∂ ϕ f ρ ∂ t = − ∇ T ρ q f (7) d f ∂ ϕ f c i ∂ t = ∇ T ϕ f D i ∇ c i − c i q f , i = C O 2   o r   H 2 S (8) q f = − d f K f μ ∇ T p − ρ g ∇ T D (9)

Where q f is the volume flow rate per unit length in the fracture, d f denotes the fracture aperture, ϕ f is the fracture porosity, K f represents the fracture permeability and ∇ T is the tangential derivatives along the interior boundary representing fractures.

COMSOL Multiphysics’ flash calculation at the thermodynamic interface was used to obtain CO₂ and H₂S concentrations at the top boundary, assuming local thermodynamic equilibrium. The Rachford-Rice equation (Rachford and Rice, 1952) is used to obtain the partitioning of different components into two phases given as: f n L = ∑ z i 1 − k i n L 1 − k i + k i = 0 (10)

Once the mole fraction of the liquid phase ( n L ) is solved, the mole fractions of components in both phases can be solved by: x i = z i n L 1 − k i + k i (11) y i = x i k i (12)

The dissolved H₂S and CO₂ compositions in the water phase are used as the top boundary conditions, as shown in Table 1. The average error between our results and the results from Li and Jiang (2020) is 1.3%, indicating that our results can be used in this study.

The domain boundaries are characterized to govern the flow and concentration behaviors. For the left, right, and bottom boundaries, impermeability is assumed, resulting in the following conditions: u x 0 , z , t = u x W , z , t = u z x , H , t = 0 (13) ∂ c ∂ x x = 0 = ∂ c ∂ x x = W = ∂ c ∂ z z = H = 0 (14)

At the top boundary, a constant pressure condition is implemented to satisfy the continuity equation: P x , 0 , t = P t o p (15)

To account for perturbations caused by CO₂ injection or formation heterogeneity, a sine wave function is introduced. This function represents the concentrations of CO₂ and H₂S at the top boundary and is defined as c x , 0 , t = 1 + A 0 ⁡ sin ⁡ ⁡ 2 π x / λ o (16)

Where A 0 is the dimensionless amplitude of the sine wave function and λ o denotes the wavelength.

Initially, the water phase contains no CO₂ or H₂S, and the pressure in the domain is equal to the pressure at the top boundary: c x , z , 0 = 0 (17) P x , z , 0 = P t o p (18)

To investigate the density-driven convection, a two-dimensional domain below the gas-saturated brine layer was built using COMSOL Multiphysics 5.6, as shown in Figure 1A. Both the domain width and height were 1 m. The formation brine in the research domain was assumed to be pure water, and no gas was initially dissolved. All domain boundaries were assumed to be closed except the top boundary, which was assumed to be a constant concentration boundary. The CO₂ and H₂S concentrations are shown in Table 1. Other assumptions included the following: First, the top boundary was assumed to be sharp, and the two-phase transition zone caused by capillary pressure was ignored (Emami-Meybodi et al., 2015). Second, the free gas mixture above the top boundary was assumed to be abundant, so the pressure change due to the dissolution of the gas mixture was ignored, and the top boundary was kept gas-saturated (Li and Jiang, 2020). Third, the domain was assumed to be isothermal (Islam et al., 2014).

To simulate the convective dissolution problem, the “Darcy’s flow” interface and the “transport of diluted species in porous media” interface were coupled. To accurately capture the diffusive behavior during the early stages, it is necessary to ensure that the grid size near the top boundary is smaller than the critical thickness of the diffusive boundary layer at the onset of convection. The critical thickness, denoted as L, can be calculated using the equation L = 3.6 c 1 μ ϕ D / g ∆ ρ k h , where c 1 is a constant factor from linear stability results (Ennis-King et al., 2005). In this study, the calculated critical thickness is 0.0076 m. To meet this requirement, the grid block size near the top boundary was set to 0.0067 m. Away from the top boundary, the grid block size gradually increased to 0.028 m in order to improve computational efficiency. The time-stepping algorithm employed in this study is the strict backward differentiation formula (BDF), with a maximum BDF order set to 2 to ensure numerical stability. To enhance computational efficiency, a segregated solver with an automatic Newton-Raphson iteration method is utilized. For the solution of the linear system, the MUMPS solver is employed, utilizing preordering algorithms to effectively permute columns and minimize fill-in. The simulation parameters used in this study are shown in Table 2. The model is validated by the published model constructed by Kim et al. (2019) on the base case without H₂S and fractures, as shown in the Supplementary Appendix.

Four measurable parameters were defined to study the dynamic behavior of convective dissolution. The first parameter is dimensionless spatial concentration, which quantifies the extent of CO₂ dissolved in water calculated by c D = ∫ 0 y ∫ 0 x c x , y d x d y C t o p ∫ 0 H ∫ 0 L d x d y (19)

This parameter is analyzed versus the dimensionless time represented by t D = k m Δ ρ g ϕ m μ H t (20)

Where k m is matrix permeability, Δ ρ is the density difference between gas-saturated water and gas-free water, g is the gravity term, ϕ m is the matrix porosity, μ is the mixture viscosity and H is the domain height.

The second parameter is the total mass flux at the top boundary ( T F t o p ), which is analyzed versus dimensionless time to describe the temporal variation of flow magnitude at the top boundary.

The third parameter is the Rayleigh number (Ra), which measures the magnitude of convective flow to diffusive flow, represented by R a = Δ ρ g k m H ϕ m μ D (21)

The fourth parameter is the Sherwood number ( S h n ) which characterizes the relative contribution of convective to diffusive flux, defined as S h n = T o t a l   f l u x   m a g n i t u d e / D i f f u s i v e   f l u x   m a g n i t u d e (22)

To study the effects of H₂S concentration on convective mixing, three H₂S concentrations (0%, 10%, and 20%) were designed. The typical features of CO₂ plume development and curve are explained using the 0% H₂S case as an example. As Figures 2, 3 show, from t D = 0 to t D = 0.21, the small fingers at the top boundary exist due to the onset of convective mixing. In this period, the T F t o p keeps increasing and reaches its peak at t D = 0.21, called the flux-growth regime. From t D = 0.21 to t D = 1.41, convective fingers start merging because of the lateral movement called the merging regime. In this regime, the T F t o p decreases from its peak. From t D = 1.41 to t D = 6.64, the complex down-swelling fingers develop and introduce a stronger counter-current around them (the white arrows around the down-swelling fingers become visible), hindering the downward movement of small fingers and making small fingers at the top boundary disappear. In this period, the merged high-concentration fingers keep moving downward, and the T F t o p stays constant, called the constant-flux regime. From t D = 6.64 to t D = 9.98, the complex down-swelling fingers hit the boundary, and the domain begins to be saturated. These down-swelling fingers then begin to move upward, taking the high-concentration solute out of the domain (the map color of the upward flow zone changes from blue to yellow), leading to the decrease in the T F t o p called the shutdown regime. At t D = 9.98, the domain is nearly saturated. The temporal regimes observed in this study, including the flux-growth regime, merging regime, constant-flux regime, and shutdown regime, align with the typical temporal regimes for convective mixing, as described by Slim (2014). This alignment serves as validation for the accuracy of the constructed model.

The density difference between gas-saturated and gas-free water, the Rayleigh number, the dimensionless CO₂ concentration at t D = 9.98, and the Sherwood number at t D = 9.98 of three H₂S concentrations are shown in Table 3. The density difference decreases with the increase in H₂S concentration. The lower density difference results in a lower convective mixing magnitude. As a result, the Rayleigh number, dimensionless CO₂ concentration, and Sherwood number all decrease as the H₂S concentration increases. The obtained results align with the conclusions drawn by Kim and Song (2017) in their linear stability analysis. They emphasized the substantial impact of impurity concentration on finger shape and growth history, with the presence of H₂S delaying the onset of convective mixing. Similarly, the results are consistent with the observations of Liu et al. (2017) in their direct numerical simulations, revealing the inhibitive influence of dissolved H₂S on convective mixing in the brine phase. The maximum Sherwood number is around 200 at t D = 9.98, indicating that the convective flow is 200 times greater than the diffusive flow, indicating that density-driven convection effectively enhances solubility trapping and significantly increases the amount of CO₂ sequestered. The lower convective mixing magnitude also delays the onset time of each flow regime shown in Figure 3. For example, the 20% H₂S case does not experience the shutdown regime at t D = 9.98, but the flux-growth regime is apparent.

This section analyzed the effects of three fracture intersection angles (30 ° , 60 ° , and 90 ° ) for three H₂S concentrations (0%, 10%, and 20%). Figure 4 compares CO₂ plume development for 10% H₂S concentration between cases with three fracture intersection angles and the case without fracture. There were two intersected fractures in the domain. One (fracture A) was horizontal and fixed. The other (fracture B) was changed to make intersection angles vary from 90 ° to 30 ° . Figures 4, 5A show that different fracture intersection angles play different roles in CO₂ plume development, resulting in an increase or decrease in the dimensionless CO₂ concentration. For the case of θ = 90 ° (10% H₂S), the presence of a highly conductive fracture B causes the plumes in the central region to merge. As a result, the merged plumes flow through fracture B and reach the bottom boundary before the plume in the left area, leading to accelerated convective mixing within the fractures. These findings are in agreement with the convective mixing analysis results for pure CO₂ presented by Kim et al. (2019). Concerning the case of θ = 60 ° (10% H₂S), the merged central plume flows in the fracture and coalesces with the left plume before hitting the bottom boundary, which slows down the accelerated convective mixing in the central region. The dimensionless CO₂ concentrations of θ = 60 ° and 90 ° are higher than that of the case without fracture, indicating that intersected fractures have a positive effect. In comparison, for the case of θ = 30 ° (10% H₂S), the fingers flowing in fracture B merge with the left plume so early that the merging effect delays the convective mixing. The dimensionless CO₂ concentration of θ = 30 ° is lower than that of the case without fracture, indicating the negative effect of intersected fractures.

The dimensionless CO₂ concentration difference between the case without fracture (10% H₂S) and the case of θ = 30 ° (10% H₂S) is prominent in the early to mid-term, and this gap is narrowed after t D = 8 (marked by the circle in Figure 5A). This is because, in the case of θ = 30 ° (10% H₂S) shown in Figure 4, the plumes in fractures flow horizontally instead of downward and merge with the left plume, which prevents the plumes from flowing across fractures in the central area in the early to mid-term. More plumes flow across fractures until convective mixing delays due to the boundary effect. Convective mixing is enhanced at this time. As shown in Figure 5A, the dimensionless CO₂ concentration of the case of θ = 60 ° (20% H₂S) is higher than that of the case of θ = 90 ° (20% H₂S), which is different from the normal trend that the dimensionless concentration of CO₂ increases with the increase of the intersection angle. This abnormal case is due to the upward flow in the fracture. Figure 5B presents the alteration of flow direction in response to changes in H₂S concentration. Notably, for the 20% H₂S concentration case, the solute exhibits an upward flow within the fracture, denoted by the black circle. Conversely, for other H₂S concentration cases, the solute flows downward. The upward flow within the fracture acts as a conduit for transporting high-concentration solute out of the domain, hindering convective mixing. Additionally, the H₂S concentration significantly impacts the onset time of flow regime in the presence of fractures. As illustrated in Figure 5A, the plume with 0% H₂S reaches the shutdown regime at t D = 8, while the plume with 20% H₂S remains in the constant-flux regime due to the decrease in density difference with increasing H₂S concentration, leading to a lower magnitude of convective mixing. The delayed onset time of each flow regime results in a decrease in the dimensionless CO₂ concentration.

In conclusion, four aspects of the influence mechanisms of intersected fractures on convective mixing are discussed. Fractures enhance convective mixing by accelerating downward solute transportation. The coalescence between merged plumes around fractures and main down-swelling plumes delays the convective mixing. The high fracture conductivity prevents plumes from flowing across fractures, hindering convective mixing. The upward flow in fractures delays convective mixing by transporting high-concentration solute out of the domain. These four mechanisms determine whether intersected fractures enhance or delay convective mixing.

The fracture conductivity ratio ( R F C D ) refers to the ratio of vertical fracture conductivity to horizontal fracture conductivity. In this section, the vertical fracture permeability was set to 10 − 8 , 10 − 9 and 10 − 10   m 2 , and the fracture conductivity ratio was 0.1, 1, and 10.

The impact of fracture conductivity ratio on CO₂ plume development in 10% H₂S concentration cases is shown in Figure 6. The solute in the vertical fracture transports relatively slowly in the case of R F C D = 0.1 (10% H₂S). Part of the solute from the vertical fracture flows into the horizontal fracture and merges with the main plumes on the left and right sides of the horizontal fracture. As the fracture conductivity ratio increases to 1, the vertical fracture transports more solute downward instead of flowing into the horizontal fracture, which makes the plume along the vertical fracture hit the bottom boundary earlier than the plume on the left side of the fracture. However, the high conductivities of vertical and horizontal fractures make the solute flow into fractures and combine with the downward plumes. This merging effect prevents plumes from flowing across the horizontal fracture (marked by the rectangle). The vertical fracture takes most of the solute in the case of R F C D = 10 (10% H₂S), making the plumes merge around the vertical fracture and flow downward. In the meantime, the strong counter-current around the vertical fracture hinders the downward movement of other fingers. As a result, the plume along the vertical fracture is large.

The comparison of dimensionless CO₂ concentrations and Sherwood numbers of three fracture conductivity ratios ( R F C D = 0.1 , 1   10 ) for all H₂S concentrations (0%, 10%, and 20%) is shown in Figures 7A, B. The dimensionless concentration before the shutdown regime increases with the rise in fracture conductivity ratio for all H₂S concentrations, which aligns with the conclusions of Rezk and Foroozesh (2019). Their study demonstrated that fractures cause fingers to deviate towards the high permeability path of the fracture, enhancing density-driven convection. The reason is that the high vertical fracture conductivity enhances the conductive mixing significantly, which can be supported by the high Sherwood number. The high magnitude of T F t o p shortens the duration of the constant-flux regime because the main plumes hit the bottom boundary earlier. Thus, the ultimate dimensionless CO₂ concentration might vary. Although an increase in fracture conductivity ratio enhances conductive mixing, it is noteworthy that the delayed onset time and reduced dimensionless CO₂ concentration resulting from an elevated H₂S concentration remain unaffected. This observation further emphasizes the crucial role of H₂S concentration in the density-driven convection problem.

The grain size sequence controlled by hydrodynamic conditions results in sedimentary formation with spatial variation in permeability. Generally, stronger hydrodynamic conditions tend to induce the deposition of coarse grains and higher permeability formations (Jin et al., 2020). The typical lithology sequences include the uniform, fining upward, fining upward with noise, fining downward, and fining downward with noise, as described in Figure 8A. In this study, the permeability contrast was set to 3. The highest permeability was 1.5 × 10 − 11   m 2 , and the lowest permeability was 0.5 × 10 − 11   m 2 . The permeability variations of the modelled lithological sequences are shown in Figure 8B. To better study the effects of lithology sequences, the average domain permeabilities were the same.

Figure 9 shows the CO₂ plume developments of different lithology sequences for 10% H2S concentration. The plume shapes are different among uniform, fining upward, and fining downward cases. The plumes of fining upward sequence do not merge until the main plumes hit the bottom boundary, which creates many thin and long fingers. This is because the downward permeability is higher than the horizontal permeability, making plumes prefer to flow downward. The lower permeability near the top boundary limits convective mixing, which causes the dimensionless CO₂ concentration of fining upward sequence to be lower than that of uniform sequence. On the contrary, the plumes of fining downward sequence start merging earlier than the plumes of uniform sequence, and strong fingers are created. The reason is that plumes tend to flow horizontally because the downward permeability is lower than the horizontal permeability. The higher permeability near the top boundary facilitates the convective mixing and increases the dimensionless CO₂ concentration. The fining upward with noise sequence accelerates the merging compared to the fining upward sequence, while the fining downward with noise sequence delays the merging compared to the fining downward sequence, given that the increased number of noise layers results in a higher permeability change gradient near the top boundary. Besides, the layered structures of CO₂ plume developments are identified in the noise sequences.

The magnitudes of dimensionless concentration, ranked from high to low, are fining downward, uniform, and fining upward lithology sequences, as shown in Figure 10A. The distinction between fining upward and upward with noise, as well as between fining downward and downward with noise, becomes more pronounced with an increase in H₂S concentration. This observation can be attributed to the interplay between density-driven convection and the effect of the noise sequence. When the H₂S concentration is low, density-driven convection is prominent, diminishing the impact of the noise sequence. Conversely, a high H₂S concentration leads to a delay in density-driven convection, amplifying the effect of the noise sequence. The curve shapes of T F t o p are different among different lithology sequences shown in Figure 10B. The curve shape of the fining upward sequence has a relatively flat peak, unlike the sharp peak in the uniform sequence. This is because the permeability increases downward, accelerating the convective mixing and extending the peak period. The curve shape of fining upward with noise sequence is basically the same as that of fining upward sequence except that the peak magnitude of fining upward with noise is higher due to the higher permeability change gradient at the top boundary. In comparison, the curve shapes of both fining downward sequence and fining downward with noise sequence have more than one peak because the formation of extra-large plumes accelerates the convective mixing. This positive effect exceeds the negative effect caused by the downward decrease in permeability.

Six large-scale cases with 20 m in width and 10 m in height were run to investigate the scale effect on density-driven convection. The relative parameters are shown in Table 4. The matrix porosity and permeability were set to 0.2 and 1 × 10 − 12   m 2 , respectively. The fracture networks in cases C and E were identical, with 16 dead-end and intersected fractures. The fracture porosity and permeability were statistically generated following the normal distribution. The mean ( μ ) and standard deviation ( σ ) of fracture porosity were 0.4 and 0.05, respectively. The mean and standard deviation of fracture permeability were 10 − 9   m 2 and 3 × 10 − 10   m 2 , respectively. The large-scale domain with a statistically generated fracture network is shown in Figure 1B. For fining upward sequence cases, the permeability increased from top ( 0.5 × 10 − 12   m 2 ) to bottom ( 1.5 × 10 − 12   m 2 ). Concerning fining downward sequence cases, the permeability decreased from top ( 1.5 × 10 − 12   m 2 ) to bottom ( 0.5 × 10 − 12   m 2 ). The average matrix permeabilities of large-scale cases were the same.

The plume developments of cases B and C are shown in Figure 11. When case B is compared to the case with 10% H₂S concentration in Figure 2, the plume developments with dimensionless time are similar, demonstrating the applicability of small-scale results to large-scale. At t D = 0.21, the plumes in both cases become apparent, and complex plumes form at t D = 1.41. The complex swelling plumes move downward and hit the bottom boundary from t D = 1.41 to t D = 6.64 for small-scale and large-scale cases. At t D = 9.98, the domain is nearly saturated, and the upward flow takes the high-concentration solute out of the domain (the map color changes from blue to yellow). In the CO₂ plume developments of case C, four main influence mechanisms of fractures on a small scale, explained in Section 3.2, can be identified on a large scale. The downward solute transportation is accelerated by fractures, which enhances convective dissolution (zone A). The coalescence between plumes around fractures and main down-swelling plumes delays convective mixing (zone B). The high fracture conductivities prevent plumes from flowing across fractures, which makes the flux magnitude below fractures low and hinders convective mixing (zone C). The flow in fractures moves upward, and high-concentration solute is transported out of the domain, delaying convective mixing (zone D).

The dimensionless CO₂ concentration and Sherwood number for six large-scale cases at t D = 9.98 are shown in Table 4. The dimensionless CO₂ concentration increases as the H₂S concentration decreases because of enhanced density-driven convection. Besides, the fining downward sequence has the highest dimensionless concentration, followed by uniform and fining upward sequences. Fractures accelerate the convective mixing in both the fining upward and uniform sequence due to the enhanced convective mixing by high conductivity fractures. The trends of Sherwood number are the same as the trends of dimensionless concentration. All trends on a large scale are the same as trends on a small scale, indicating the applicability of small-scale results to a large scale.