Enhanced steady-state dissolution flux in reactive convective dissolution

TL;DR

This paper investigates how chemical reactions influence the development and intensity of convective dissolution in stratified fluids, revealing that reactions can enhance steady-state flux and convection, with implications for CO2 sequestration.

Contribution

It provides a numerical analysis of reactive convective dissolution dynamics, highlighting how reactions modify flux and convection intensity depending on chemical parameters.

Findings

Reactions can accelerate or slow down convection.

Chemical reactions enhance steady-state dissolution flux.

Reactions can induce more intense convection than non-reactive cases.

Abstract

Chemical reactions can accelerate, slow down or even be at the very origin of the development of dissolution-driven convection in partially miscible stratifications, when they impact the density profile in the host fluid phase. We numerically analyze the dynamics of this reactive convective dissolution in the fully developed non-linear regime for a phase A dissolving into a host layer containing a dissolved reactant B. We show that for a general A+B$\rightarrow$C reaction in solution, the dynamics vary with the Rayleigh numbers of the chemical species, i.e. with the nature of the chemicals in the host phase. Depending on whether the reaction slows down, accelerates or is at the origin of the development of convection, the spatial distributions of species A, B or C, the dissolution flux and the reaction rate are different. We show that chemical reactions enhance the steady-state flux as they consume A and can induce more intense convection than in the absence of reactions. This result is important in the context of CO$_2$ geological sequestration where quantifying the storage rate of CO$_2$ dissolving into the host oil or aqueous phase is crucial to assess the efficiency and the safety of the project.

Figures

Tables

Table 1. Table 1: Specific reactive cases discussed in 3 : R1, R2 and R3. In all cases β = 1 𝛽 1 \beta=1 .

	-1	1
1	R1 less unstable reactive	R2 more unstable reactive
-1	(not discussed) stable	R3 unstable due to reaction

Equations

$$A$$

$$t$$

$$\frac{\partial A}{\partial t}+(\mathbf{u}\cdot \nabla )A$$

$$\frac{\partial B}{\partial t}+(\mathbf{u}\cdot \nabla )B$$

$$\frac{\partial C}{\partial t}+(\mathbf{u}\cdot \nabla )C$$

$$u_z(z=0)=0;$$

$$A(z=0)=1;$$

$${\frac{\partial B}{\partial z}|}_{z=0}=0;$$

$${\frac{\partial C}{\partial z}|}_{z=0}=0;$$

$$A(x,z=0,t=0)=1+ϵ\cdot \mathrm{rand}(x);A(x,z>0,t=0)=0,$$

$$B(x,z,t=0)=\beta =B_0/A_0,$$

$$C(x,z,t=0)=0,$$

$$\tilde{\rho }={\rho }_0(1+{\alpha }_A\tilde{A}+{\alpha }_B\tilde{B}+{\alpha }_C\tilde{C}),$$

$$\frac{\tilde{\rho }-{\rho }_0}{{\rho }_c}=R_{A}A+R_{B}B+R_{C}C,$$

$$R_i=\frac{{\alpha }_i{A}_{0}g\kappa l_c}{\varphi \nu D_A}=\frac{{\alpha }_i{A}_{0}g\kappa }{\varphi \nu \sqrt{D_{A}q{A}_0}},$$

$$B=\beta -C.$$

$$\rho =\frac{\tilde{\rho }-{\rho }_0}{{\rho }_c}-R_B\beta ,$$

$$\rho =R_{A}A+\Delta R_{CB}C,$$

$$\nabla p$$

$${\nabla }^{2}p$$

$$\bar{r}(z,t)=1/L{\int }_0^{L}A(x,z,t)B(x,z,t)\mathrm{d}x.$$

$$\bar{f}(z,t)=1/L{\int }_0^{L}f(x,z,t)\mathrm{d}x,\text{with f=A,B,C,}$$

$$z_m(t)=2\mathrm{erfinv}(1-s)\sqrt{t},$$

$$z_m(t)=2\mathrm{erfinv}(1-s/2)\sqrt{t},$$

$$\bar{\lambda }(t)={\left(\frac{{\int }_{1/L}^{1/(2dx)}\nu |\mathcal{F}(\bar{\rho }){|}^2\mathrm{d}\nu }{{\int }_{1/L}^{1/(2dx)}|\mathcal{F}(\bar{\rho }){|}^2\mathrm{d}\nu }\right)}^{-1},$$

$$z_f=2{\eta }_f\sqrt{t},$$

$$z_f=\frac{{\int }_0^{H}z\bar{r}(z)\mathrm{d}z}{{\int }_0^H\bar{r}(z)\mathrm{d}z},$$

$$w_f=2{\left(2\ln (1/\xi )\frac{{\int }_0^H(z-z_f{)}^2\bar{r}(z)\mathrm{d}z}{{\int }_0^H\bar{r}(z)\mathrm{d}z}\right)}^{1/2}.$$

$$⟨c_i⟩=\frac{1}{V}{\int }_{\Omega }{c}_i\mathrm{d}V.$$

$$\frac{\partial ⟨A⟩}{\partial t}=\frac{J}{H},$$

$$\frac{\partial ⟨A⟩}{\partial t}$$

$$\frac{\partial ⟨C⟩}{\partial t}$$

$$\frac{\partial ⟨A+C⟩}{\partial t}$$

$$J=-\frac{1}{L}{\int }_0^L{\frac{\partial A}{\partial z}|}_{z=0}\mathrm{d}x.$$

$$J_D=\frac{1}{\sqrt{\pi t}},$$

$$J_{RD}=\frac{1+\beta }{\sqrt{\pi t}},$$

Full text

[TABLE]

00footnotetext: ∗* Corresponding author; E-mail: adewit@ulb.ac.be00footnotetext: a Université libre de Bruxelles (ULB), Faculté des Sciences, Nonlinear Physical Chemistry Unit, CP231, 1050 Brussels, Belgium.00footnotetext: b Université libre de Bruxelles (ULB), Faculté des Sciences, Fluid and Plasmas Dynamics Unit, CP231, 1050 Brussels, Belgium.*

1 Introduction

Dissolution-driven convection can develop in partially miscible stratifications when the dissolution of a phase A with a finite solubility into a host fluid phase creates an unstable density stratification. This can happen for instance when the phase A, dissolving from above, increases the density of the host phase, thereby forming a layer of denser fluid rich in A on top of less dense fluid. Studying such convective dissolution can help improve the safety of nuclear reactors 1, optimize industrial production of chemicals 2, 3 or understand the physicochemical processes at hand during CO2 geological sequestration 4, 5. The temporal evolution of the dissolution-driven convective dynamics has been characterized in detail and it has been shown that the dissolution flux reaches a steady-state value before shutdown 6, 7, 8, 9, 5. Understanding the impact of chemical reactions on such dynamics has recently gained interest because of the potential effect of geochemistry on the efficiency of CO2 geological sequestration 5, 10. Reactions can indeed affect convection because they modify solute concentrations affecting fluid properties such as the density of the solution, leading to different possible scenarios for the development of convection 11, 12.

When the dissolving species A reacts with a solute B initially present in the host solution following a second-order scheme A+B $\rightarrow$ C, the reaction can slow down or accelerate the development of convection compared to the non-reactive case depending on the type of density profile building up in the host phase, as shown by both experimental and numerical approaches 13, 14, 15, 16, 17, 18, 2, 3. If C is less dense than B, a density profile with a minimum is observed and the instability develops more slowly than in the non-reactive case. On the contrary, if C is sufficiently denser than B, the instability develops faster and the density profile is monotonic like its non-reactive counterpart. A simpler reaction where B is solid and no C is produced reduces the growth rate of convection, because A at the origin of the instability is consumed by the reaction 19, 20, 21, 22, 23, 24, 25, 26, 27, 28. In contrast, when A decreases the density of the solution upon dissolution, the non-reactive case is buoyantly stable and reactions can be at the origin of a density profile unstable with regard to convection due to the creation of a maximum of density 14, 18.

The different possible convective dynamics in the presence of reaction have been classified according to whether the reaction slows down or accelerates the development of dissolution-driven convection 13, 14, 18, 16. We here broaden the scope of this classification by analyzing the effects of reaction on other aspects of the convective dynamics in the fully developed non-linear regime. First, it is crucial to understand whether the different successive regimes of the convective dynamics identified in the non-reactive case can also be observed in reactive cases. Second, for potential applications it is of interest to characterize the impact of reactions on the evolution of the flux of A dissolving into the host phase, in particular on its steady-state value already quantified in the non-reactive case 6, 7, 29, 8, 9, 30, 31, and on the global reaction rate in the host phase. Finally, we aim to characterize the effects of convection on the dynamics of the reaction front, which also remain poorly understood.

We address these issues by theoretically studying the effects of an A+B$\rightarrow$C reaction on the non-linear dynamics during dissolution-driven convection in a host fluid phase occupying a porous medium, relevant to the context of CO2 sequestration. Similarly to previous studies 13, 14, 15, 32, 17, we consider equal diffusivities of the three chemical species in order to focus solely on solutal effects. By varying the relevant control parameters, we study their effect on the properties of the reaction-diffusion-convection dynamics: the evolution of the fingering pattern, the dynamics of the reaction front, the dissolution flux and the global reaction rate affecting the storage of A into the host fluid are all quantified. We qualitatively describe the dynamics observed in the unstable non-reactive case and in three typical reactive cases: less unstable or more unstable than the unstable non-reactive counterpart, as well as destabilization by reaction of a stable non-reactive counterpart. Further, we derive simple scalings to predict the steady-state reaction rate and dissolution flux of A at the interface as a function of the Rayleigh numbers. Such scalings could be used to predict the temporal evolution of the quantity of A in the host phase.

2 Model

We consider an isothermal, isotropic and homogeneous vertical system, in which two partially miscible phases are placed in contact along a horizontal interface in a porous medium 13, 14. The gravity field $\bm{g}$ points downwards, along the vertical $\tilde{z}$ axis perpendicular to the horizontal $\tilde{x}$ axis. Phase A dissolves into the other lower fluid phase, called “host” phase as shown in Fig. 1. The concentration of A at the interface is considered to remain constant over time, and to be equal to its solubility $A_{0}$ in the host phase, following the assumption of local chemical equilibrium.

The host phase contains a reactant B dissolved with an initial concentration $B_{0}$. Species A reacts with solute B to produce another solute C, following a second-order A+B$\rightarrow$C reaction. All three species can thus contribute to changes in density. The concentrations of B and C are assumed small enough and therefore do not significantly affect the solubility $A_{0}$ of A into the host phase.

The interface is considered to be permeable to species A but impermeable to the solvent of the host phase and to solutes B and C. To focus on the effects of the reaction on the dynamics, we assume that the interface remains in the course of time at the same position $\tilde{z}=0$ and we study the dynamics in the host phase only. The host phase extends from $\tilde{x}=0$ to $\tilde{L}$ in the horizontal direction and from the interface at $\tilde{z}=0$ to $\tilde{z}=\tilde{H}$ in the vertical direction.

An equation for the fluid flow velocity $\bm{\tilde{u}}=(\tilde{u}_{x},\tilde{u}_{z})$ of an incompressible flow is coupled to the reaction-diffusion-convection (RDC) equations for solute concentrations via an equation of state for the density depending linearly on the concentrations13, 14. The concentrations, time, spatial coordinates and velocity are normalized using 13, 14:

[TABLE]

where tildes denote dimensional variables. We nondimensionalize the solute concentrations $\tilde{A}$, $\tilde{B}$ and $\tilde{C}$ in Eq. (1) with the solubility $A_{0}$ of phase A. In Eq. (2), we use the chemical time scale $t_{c}=1/(qA_{0})$ with $q$ the kinetic constant of the reaction A+B$\rightarrow$C, the RD length scale $l_{c}=\sqrt{D_{A}t_{c}}=\sqrt{D_{A}/(qA_{0})}$ with $D_{A}$ the diffusion coefficient of A and the velocity scale $u_{c}=\phi l_{c}/t_{c}=\phi\sqrt{D_{A}qA_{0}}$ with $\phi$ the porosity of the medium at hand.

The dimensionless RDC equations then read:

[TABLE]

where $\delta_{B}=D_{B}/D_{A}$ and $\delta_{C}=D_{C}/D_{A}$, with $D_{B}$ and $D_{C}$ the diffusion coefficients of species B and C. Using the notations $H=\tilde{H}/l_{c}$ and $L=\tilde{L}/l_{c}$ for the dimensionless height and width of the host phase, we solve Eqs. (3)-(5) with the boundary conditions (see Fig. 1)

[TABLE]

and the initial conditions

[TABLE]

where $\beta=B_{0}/A_{0}$ is the ratio between the initial concentration $B_{0}$ of reactant B and the solubility $A_{0}$ of A in the host phase. When $\beta=0$, the non-reactive case is recovered; when $\beta\to\infty$, reactant B is in large excess with regard to A so that the reaction can be considered first-order at early times 21, 22, 23, 24, 25, 26, 27, 32. Equation (10) expresses that perturbations are introduced in the initial concentration of A at the interface in order to trigger the instability (see e.g. Refs. 33, 34 for a discussion of the possible types of perturbations). $\epsilon\ll 1$ is the amplitude of the perturbation, here chosen as $10^{-3}$, and rand$(x)$ is its modulation, function of the horizontal coordinate $x$ and varying randomly between -1 and 1 (“white noise”).

The set of equations (3)-(5) are closed using an equation for the fluid flow velocity of an incompressible flow. To that end, we assume a linear state equation for the dimensional density $\tilde{\rho}$ of the solution as a function of the solute concentrations:

[TABLE]

where $\alpha_{i}=\frac{1}{\rho_{0}}\frac{\partial\rho}{\partial c_{i}}$ is the solutal expansion coefficient of species i. A dimensionless density can be computed as

[TABLE]

where $\rho_{c}=\phi\mu D/(g\kappa l_{c})$ is the density scale and the Rayleigh numbers $R_{i}$ (i=A,B,C) quantify the contribution of species i to the dimensionless density of the solution, constructed with the RD length scale (Eq. (2)):

[TABLE]

with $\nu=\mu/\rho_{0}$ the kinematic viscosity of the solvent. The expression (14) is appropriate in the general case where the species have different diffusivities. For $\delta_{i}\neq 1$, the problem is thus dependent on six parameters: $\delta_{B}$, $\delta_{C}$, $R_{A}$, $R_{B}$, $R_{C}$ and $\beta$.

We further assume here that all species A, B and C have the same diffusion coefficient so that $\delta_{B}=\delta_{C}=1$. We can then add Eqs. (4) and (5), taking into account the corresponding boundary (Eqs. (8)-(9)) and initial (Eqs. (11)-(12)) conditions, to obtain the conservation relation

[TABLE]

Thanks to this conservation relation, we can reduce the number of dimensionless parameters further by defining the dimensionless density $\rho$ as

[TABLE]

i.e.

[TABLE]

where $\Delta R_{CB}=R_{C}-R_{B}$ represents the difference between the contributions to density of product C and reactant B. For equal diffusivities $\delta_{B}=\delta_{C}=1$, the definition 18 of $\rho$ explicitly highlights that the system is characterized by only three parameters, here chosen as111Note that our previous results13, 14 remain the same with this new formulation 18 for $\rho$ because they were performed for $R_{B}=0$, so that the parameter $R_{C}$ that we varied is strictly equivalent to $\Delta R_{CB}$ here. $R_{A}$, $\Delta R_{CB}$ and $\beta$. Darcy’s equation and the incompressibility condition expressed as Poisson equation then read in dimensionless form:

[TABLE]

with $p$ the dimensionless pressure.

We numerically solve Eqs. (3)-(5) with $\delta_{B}=\delta_{C}=1$, Eqs. (18)-(20) with Eqs. (6)-(12) on a computational domain of width $L=3072$ and height $H=2048$ using the YALES2 software 35, more specifically the DARCY_SOLVER module 36. This software is based on the finite volume method 37. We use an explicit method called TFV4A or TRK4 (two-step Runge-Kutta with a fourth-order spatial discretization) 38. The dynamics depend on the random noise added to the initial condition in Eq. (10). Therefore, for each value of the set of parameters ($R_{A},\Delta R_{CB},\beta$), we average the results over 15 realizations to obtain robust results. Increasing that number of realizations above 15 does not impact the averages and standard deviations of the results significantly (below 5%). The uncertainty linked to the different possible noises is quantified as the 95 % confidence interval for two-sided critical regions. In addition, we have checked that these results averaged over 15 realizations were robust with regard to refinement of the iterative convergence tolerance (here 10*-10*) for solving Poisson’s equation (20) with HYPRE, mesh size ($\Delta x=\Delta z=4$) and time step ($\Delta t=0.5$). With these values, the iterative convergence errors and discretization errors on the results were smaller than 5%.

We perform a parametric study of the non-linear dynamics as a function of $R_{A}$ and $\Delta R_{CB}$ while keeping $\beta=1$. Indeed it is already known that increasing $\beta$, i.e. amplifying the amount of dissolved reactant B with regard to the solubility of A, amplifies the effect of reaction on the development of dissolution-driven convection. If the instability develops faster than in the non-reactive case, a larger $\beta$ accelerates even more that development; conversely if chemistry slows down the growth of convection, increasing $\beta$ decreases the growth rate of the instability 13, 14, 15, 16, 17. Our objective here is to analyze the effect of changing the nature of reactants on dissolution-driven convection. We therefore analyze two main classes of dissolving species A, taking $R_{A}$=+1 (-1) representing a component increasing (respectively decreasing) the density of the host phase upon dissolution. Fixing thus $|R_{A}|=1$ and only varying the sign of $R_{A}$ affects the type of density profile above the reaction front 12, 14: increasing downwards if $R_{A}=-1$ (stable non-reactive counterpart) or decreasing downwards if $R_{A}=1$ (unstable non-reactive counterpart). For each case, we vary $\Delta R_{CB}$, which corresponds to scanning various possible reactant B and product C pairs. To study the effect of the composition and, in particular, of solutal effects on the convective dynamics, we analyze the fingering and reaction zone evolution, the dissolution flux, and the volume-averaged concentrations as a function of $\Delta R_{CB}$ varying between -1 and +1.

3 Dissolution-driven convective dynamics

We start by qualitatively describing the dissolution-driven convective dynamics in a few specific cases without reaction (Section 3.1) or with reaction (Section 3.2). We then compare the fingering dynamics using space-time plots (Section 3.3).

3.1 Non-reactive cases

The non-reactive case can be either unstable ($R_{A}>0$) or stable ($R_{A}<0$) with regard to dissolution-driven convection.

3.1.1 NR case: non-reactive unstable ($R_{A}=1$)

The non-reactive (NR) unstable case for which $R_{A}=1$ has been well characterized in the literature 6, 7, 8, 9, 5. If species A increases the density of the host phase ($R_{A}>0$), the dissolution of A into a lower host phase progressively creates a buoyantly unstable density stratification. The dynamics in this unstable non-reactive case are illustrated for a specific realization in Fig. 2 showing the density field at different times. The dynamics can be divided in different successive regimes 8, 9. Initially, the miscible contact zone between the denser zone rich in A below the interface and the less dense bulk solution below it is flat (Fig. 2a). This zone deforms gradually once fingers of the denser fluid begin to sink into the lower part of the host phase (Fig. 2b). At the beginning these fingers do not interact significantly with each other and a well defined wavelength is observed: soon after the onset of the instability, 22 fingers can be observed on the total width of 3072 (see Fig. 2b), which corresponds to a wavelength of $\approx$ 140, in agreement with the results of other non-linear simulations 9, 7. After some time, merging becomes a dominant process and the number of fingers decreases dramatically (Figs. 2c-2d). After this merging regime, the number of fingers only decreases slightly (Figs. 2e-2f). In the reinitiation regime 39, 9, small new fingers, called protoplumes, develop from the boundary layer and then join older fingers.

3.1.2 NR2 case: non-reactive stable ($R_{A}=-1$)

If the density decreases upon dissolution of A ($R_{A}\leq 0$), the stratification is buoyantly stable in the absence of reactions as the growing boundary layer rich in A is less dense than the host solvent. Species A then invades the host phase by diffusion only, and no convection develops (not shown here).

3.2 Reactive cases

We now analyze the effect of chemical reactions on the two non-reactive cases presented above. To do so, we describe three specific reactive cases as shown in 1.

3.2.1 R1 case: less unstable reactive system ($R_{A}=1$, $\Delta R_{CB}=-1$)

When the solution of C is less dense than that of B ($\Delta R_{CB}<0$), the chemical reaction is expected to slow down the development of fingering because a minimum develops in the density profile at the location of the reaction front 15, 16, 13, 14. The temporal dynamics of the density field shown in Fig. 3 indeed illustrate that fingers develop more slowly than in the NR case. The fingering dynamics can be described by the same successive regimes as in the NR case: fingers have not developed yet (Fig. 3a), they grow without interactions between them (Fig. 3b) and they merge several times (Figs. 3c-3f). By contrast to the non-reactive counterpart, protoplumes already form in the merging regime (Fig. 3f). The initial number of fingers (see Fig. 3b) is larger than in case NR (Fig. 2b), i.e. we have here roughly 26 fingers.

Figures 4 shows the concentration fields responsible for the density field in Fig. 3f. The dissolving species A is seen to be at the origin of the fingering pattern (Fig. 4a). There are no fingers in the concentration fields of reactant B and product C although the contact line between the reacted (where A and C are present) and unreacted (with mostly B) zones is deformed by the finger tips of A (Figs. 4b-4c). The reaction occurs mostly in a thin zone localized at the tip of the fingers (Fig. 4d), which corresponds to the position of the minimum of density as shown by the comparison of the reaction rate map (Fig. 4d) with the density field (Fig. 3f).

To analyze how the reaction zone evolves in time, we compute the horizontally-averaged reaction rate profile as

[TABLE]

The reaction profiles shown in Fig. 5a illustrate that the localized reaction zone enlarges and moves downwards in time. In parallel, the value of the maximum reaction rate decreases progressively. The reaction zone initially consists of a single peak that progressively deforms (see for instance $t=20000$ where two local maxima are visible). This can be explained by the intense coalescence occurring around that time: fingers have different lengths as merging fingers are longer than their neighbours (Fig. 3d). The reaction zone is initially symmetric, but this symmetry is lost in time as its tail enlarges more than its head, where fingers are arrested by the minimum of density.

All these characteristics of the reaction zone can be linked to those of the horizontally-averaged concentration profiles, computed as

[TABLE]

and illustrated at different times in Fig. 5b-c. At early times < 8000, the $\bar{A}$ profile (Fig. 5b) looks like its RD counterpart 14, 13, 12 as convection is not large enough to significantly affect the concentration profiles, while later $\bar{A}$ is deformed due to the apparition of fingering. At $t=28000$, we observe a plateau between $z\approx 80$ and 420 in which $\bar{A}\approx 0.3$. On a given distance, $\bar{A}$ remains almost constant as fingers contain the same quantity of A, but are thinner near the interface and more spread out just above the minimum of density, which acts as a barrier that prevents fingers from progressing further downwards in the solution (Fig. 4a). The plateau of $\bar{A}$ ends where $\bar{B}=\beta-\bar{C}$ starts to increase (i.e. $\bar{C}$ starts to decrease, see Fig. 5c) as species A is consumed by the reaction and $\bar{r}$ rises (Fig. 5a). Lower in the solution, $\bar{r}$ reaches a maximum before decreasing as A has not diffused far enough and has been depleted by the reaction with B (Fig. 5a). Similarly to the A profile, the C profiles initially look like their RD counterparts 13, 14, 12, i.e. error function-like curves which decrease from their maximum value $\beta=1$ at the reaction front to zero in the bulk of the solution (Fig. 5c). However, no bumps appear in $\bar{C}$ because the fingering pattern is mainly due to the denser dissolving species A (Fig. 4). Analyzing the density profile $\bar{\rho}(z,t)$ plotted in Fig. 5d shows that the width of the minimum of density enlarges progressively, so that it transitions in time from a strict local minimum to a zone where the density is constant, corresponding to the plateau value in $\bar{A}$ at $t=28000$ (Fig. 5a).

3.2.2 R2 case: more unstable reactive system ($R_{A}=1$, $\Delta R_{CB}=1$)

When the contribution of the product C to the density is sufficiently larger than that of the dissolved reactant B, reactions destabilize the system even more, i.e. increase the characteristic growth rate in the linear regime 15, 13, 14, 16. Similarly, fingers develop more quickly and elongate more rapidly than in the NR case (Fig. 6). The fingering dynamics in case R2 are characterized by the same successive regimes as in case NR: no fingering (Fig. 6a), linear finger growth (Fig. 6b), merging (Figs. 6c-6d), and reinitiation (Figs. 6e-6f). The number of fingers is initially larger than for NR case (compare Fig. 6b with Fig. 2b), in agreement with previous experimental studies 17, 40, 15, 16, 13 and theoretical studies 15, 14, 16.

We now analyze in Fig. 7 the concentration fields at a given time to understand how they combine to form the fingering pattern observed in the density field (Fig. 6d). The dissolving species A does not penetrate far in the host solution (Fig. 7a). The reactant A is indeed readily consumed as soon as it enters the host phase, which limits its progression into the host solution, in agreement with recent results 17, 32. The reaction must be fed by B diffusing to the upper part of the solution, close to the interface, where it reacts with dissolved A (Fig. 7b). The fingering pattern observed in the density field (Fig. 6d) is mainly due to the denser product C produced close to the interface, sinking into the lower part of the host solution (Fig. 7c). The convection mechanism is thus essentially the same as in the non-reactive counterpart (NR case), except that the contribution of the denser product C adds to that of the dissolving species A. This explains why the same types of regimes are observed in both cases. The reaction rate $AB$ is the largest along the side of the birth zone of the density fingers and the contact line between the boundary layer rich in A and the bulk solution (Fig. 7d).

3.2.3 R3 case: unstable due to reaction ($R_{A}=-1$, $\Delta R_{CB}=1$)

When $R_{A}<0$, the non-reactive case is stable and A+B $\rightarrow$ C reactions can be at the origin of buoyancy-driven convection as soon as $\Delta R_{CB}>0$ because a maximum, corresponding to a locally unstable stratification, then develops in the density profile 14, 12. Therefore we expect a stable boundary layer just below the interface followed further away at a given distance from the interface by a locally unstable zone generating fingers sinking down. Our numerical results confirm this prediction as shown in Fig. 8. Fingers develop from the maximum of density located below the interface and their base has a characteristic shape not observed in the other cases. Apart from that, the dynamics are similar to those observed in case R2 (Fig. 6): no fingering (Fig. 8a), linear finger growth (Fig. 8b), merging (Figs. 8c-8d), and reinitiation (Figs. 8e-8f). The reinitiation mechanism can be observed in Fig. 8f as “pulses” in the fingers, corresponding to merging with protoplumes.

Like in R2 case, the fingering pattern (Fig. 8d) is also mostly formed by solute C (Fig. 9c) sinking into the less dense solution of B (Fig. 9b) while most of the dissolving species A is consumed as soon as it enters the solution (Figs. 9a-9d). The remaining A on top decreases the density of the solution and is thus at the origin of the specific shape of the fingers’ base. Except that specific fingers’ shape, the dynamics is thus qualitatively similar but slower than that in R2 case, probably due to the stable boundary layer between the interface and the maximum of density.

3.3 Comparison of the space-time plots

The space-time plots shown in Fig. 10 summarize the differences between the dynamics in the NR case and those in the three specific reactive cases R1, R2 and R3 discussed above. These pictures are constructed by plotting the density along a horizontal line at $z=64$, except for case R3 where we had to take a line below at $z=128$ because of the stable boundary layer close to the interface. The dynamics in case R1 (Fig. 10b) are slower while those in case R2 (Fig. 10c) are faster than in case NR (Fig. 10a) because of the different density profiles in the host phase. In particular, fingers appear earlier and new protoplumes are generated more frequently in case R2 than in case NR, and vice versa for case R1. In case R3, fingering occurs (Fig. 10d) because of the effect of the reaction on the density profile. The successive regimes describing the dynamics are similar for cases NR, R2 and R3, while in case R1, some reinitiation already starts during the merging regime. For all reactive cases R1, R2 and R3, the wavelength of the fingering pattern when fingers first become visible is smaller than for case NR, in agreement with theoretical and experimental predictions 13, 14, 15, 16, 17, 40.

4 Effect of reaction-diffusion-convection interplay on the spatio-temporal dynamics

Motivated by the differences between the dynamics of cases NR, R1, R2 and R3, we further analyze the spatio-temporal dynamics of fingering and of reaction when changing $\Delta R_{CB}$ between -1 and 1. Our aim is to quantify how the reaction affects the convective dynamics and conversely how convection impacts the dynamics of the reaction zone. We have previously classified the effects of reaction on convection by evaluating a characteristic growth rate in the linear regime 13, 14. Below a given value $\Delta_{R}$ for $\Delta R_{CB}$, this growth rate is smaller than its non-reactive counterpart, and conversely. We now revisit this classification by examining various aspects of the convective dynamics in the fully developed non-linear regime before shutdown: firstly in Fig. 4 the spatio-temporal dynamics, i.e. the evolution of the fingering pattern (elongation, wavelength) and of the reaction zone, and secondly in Fig. 5 the storage properties, i.e. at what rate and under which form A is stored into the host fluid phase. Note that the trends observed when varying $\Delta R_{CB}$ for a given $R_{A}$, i.e. when changing the reactant B – product C pair for a given dissolving species A, are similar for $R_{A}>0$ and $R_{A}<0$. These trends will therefore be illustrated for $R_{A}=1$ only as the same discussion can be repeated for $R_{A}=-1$.

4.1 Mixing length and velocity

We first analyze the impact of changing $\Delta R_{CB}$ on the evolution of the mixing length and on the finger velocity. We define the mixing length $z_{m}$ as the most advanced position along $z$ where $A+C>s$, with $s$ a small arbitrary threshold here chosen as $0.01$. This position evolves dynamically in time as A dissolves into the host solution and its reaction with solute B produces C. This definition of the mixing length represents the extension of the zone containing stored A in the form of either dissolved unreacted A or product C (i.e. reacted A).

We can derive analytical expressions for the mixing length valid before convection sets in. In the non-reactive case, inserting $A(z_{m})=s$ into the diffusive concentration profile 41, 42 ${\rm erfc}(z/(2\sqrt{t}))$ of A gives

[TABLE]

$\approx 3.64\sqrt{t}$ for $s=0.01$. In the reactive case, we introduce $C(z_{m})=s$ into the RD concentration profile 14 $2{\rm erfc}(z/(2\sqrt{t}))$ of C valid for $\beta=1$ below the reaction front to get

[TABLE]

$\approx 3.97\sqrt{t}$ for $s=0.01$. Equations (23)-(24) show that the mixing length in the diffusive regime increases in time as $\sqrt{t}$, and increases faster in the presence of a reaction.

The mixing lengths computed from the numerical simulations are shown in Fig. 11. One curve represents the average over 15 realizations. The 95% confidence interval shown as lighter areas around the curves represents the variability due to the random noise on the initial condition. We see that the variability between realizations is amplified when $\Delta R_{CB}$ increases, because of more intense convection.

Both expressions (24) and (23) (shown in Fig. 11 as dotted and dashed curves respectively) are valid as long as diffusion remains the dominant transport process. After a certain transition time noted $t_{NL}$, convection becomes important, fingers start to move faster and $z_{m}$ starts departing from the $\sqrt{t}$ curve. $z_{m}$ then increases approximately proportional to $t$, although fingers might progressively slow down in some cases (see e.g. $\Delta R_{CB}=0.5$ after $t=7500$ in Fig. 11), or exhibit two successive different velocities. $t_{NL}$ represents the time when non-linearities significantly affect the vertical elongation of the fingers. More precisely, we evaluate $t_{NL}$ as the time when the relative difference between $z_{m}$ and the diffusive prediction (Eq. (23) for non-reactive cases or Eq. (24) for reactive cases) becomes larger than 5%. This time $t_{NL}$ decreases when $\Delta R_{CB}$ increases (Figs. 11, 12a), which is coherent with the predictions of the linear stability analysis that the destabilizing effect of chemistry increases when $\Delta R_{CB}$ increases 14. For $\Delta R_{CB}\leq 0.1$, $z_{m}$ deviates from the diffusive curve later than in the non-reactive case. This value is of the same order as the critical value 0.32 predicted by linear stability analyses 13, 14. The difference might arise because $t_{NL}$ is measured here on the basis of fingering dynamics rather than on the basis of the perturbation with regard to the base state.

We therefore also define a time $t_{0}$ for the onset of the instability on the basis of the perturbation in velocity, computed as $U^{2}(t)=\int_{0}^{H}\int_{0}^{L}u_{x}^{2}(x,z,t)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}z+\int_{0}^{H}\int_{0}^{L}u_{z}^{2}(x,z,t)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}z$. To highlight the dynamics at early times, we plot $U^{2}(t)$ in a log-log graph (Fig. 13). For any value of $\Delta R_{CB}$, the perturbation initially decreases until a given onset time $t_{0}$ when it reaches its minimum. Like $t_{NL}$, this onset time $t_{0}$ increases when $\Delta R_{CB}$ is increased (Fig. 12b). In other words, both linear growth and non-linear regimes start earlier if $\Delta R_{CB}$ is larger. The critical value $\Delta R_{CB}=0.2$ above which the reaction makes the system more unstable is slightly larger than the one calculated with $t_{NL}$, but still smaller than the prediction of the linear stability analysis 13, 14.

In addition, we classify the effects of reaction on convection at later times, based on the time-averaged finger velocity $\dot{z}_{m}$ computed as the least-squares fitted slope of $z_{m}(t)$ between $t_{NL}$ and the time when the fingers arrive at 95% of the depth of the host phase. $\dot{z}_{m}$ increases when $\Delta R_{CB}$ increases (Fig. 12c), so that as above we can define a critical value of $\Delta R_{CB}$ above which fingers advance faster than their non-reactive counterparts. This critical value of 0.5 is larger than that evaluated on the basis of $t_{NL}$ and $t_{0}$. Indeed for $\Delta R_{CB}=0.5$, although fingers start to accelerate earlier than their non-reactive counterparts (smaller $t_{NL}$), they move more slowly at later times and arrive at the bottom at the same time as in the non-reactive case (Fig. 11). This can be explained by non-linearities and interactions between fingers that can slow down the vertical progression of their tips in the host solution. In addition, for all measurements $t_{NL}$, $t_{0}$ and $\dot{z_{m}}$, convection starts earlier and fingers progress more slowly when $R_{A}=-1$ than for $R_{A}=1$, probably because of the buoyantly stable zone between the interface and the reaction front.

In summary, classifying the effects of reactions on convection can be based on different criteria. This classification depends on whether we are interested in the onset of the convective instability, in the time when convection becomes visible or in the average progression of the fingers in the host fluid. For intermediate values of $\Delta R_{CB}$ between 0.1 and 0.5, convection can indeed start earlier but fingers still progress more slowly than in the non-reactive case. This result highlights that the criterium chosen for the classification depends on what is required for the application: earlier developing or more intense convection.

4.2 Wavelength of pattern

The power-averaged mean wavelength $\bar{\lambda}$ of the fingering pattern is computed as 43, 44

[TABLE]

where $\nu=1/\lambda$ is the number of wavelengths per unit distance and $\mathcal{F}(\bar{\rho})$ is the Fourier transform of the vertically averaged density profile $\bar{\rho}(x,t)$ evaluated as $\bar{\rho}(x,t)=\frac{1}{H}\int_{0}^{H}\!\rho(x,z,t)\mathop{}\!\mathrm{d}z$. In the following figures, the results are shown until the start of the shutdown regime occurring after fingers have touched the bottom of the solution. The wavelength $\bar{\lambda}$ typically increases over time as shown in Fig. 14, i.e. the number of fingers decreases as fingers become wider or merge with each other. This increase follows a $\sqrt{t}$ trend in the first regimes of the convective dynamics, i.e. diffusive, linear growth and flux growth, when the roots of the fingers remain mostly immobile. At later times corresponding to the merging regime, $\bar{\lambda}$ increases approximately linearly with time. During this linear increase, two different slopes can be distinguished: the first one is larger, expressing intense merging, while the second one is smaller. For example in case NR, between $t\approx 7500$ and 10000 the wavelength increases from 200 to 450 over this time period of 2500, while after $t\approx 10000$, the wavelength increases only up to 750 over a period of 15000 (see Fig. 10). The slope of $\bar{\lambda}(t)$ increases with $\Delta R_{CB}$, meaning that fingers merge faster. In particular, the intense merging can occur faster (approximately when $\Delta R_{CB}\geq 0.5$) or more slowly (for $\Delta R_{CB}\leq 0.2$) than in NR case.

By comparing Figs. 14 and 11, we note that the vertical and horizontal dynamics are linked. In the merging regime, fingers advance more slowly but merge faster. When two fingers merge, the resulting finger is denser as solute coming from two fingers has accumulated in only one finger, a more confined region of space. Therefore, during the next stage, the resulting denser fingers sink faster. They also merge more slowly, probably due to their amplified velocity. However, for $\Delta R_{CB}\leq 0.2$, fingers do not accelerate after intense merging, probably because the consequent increase of the average finger weight is not large enough to significantly affect the vertical finger velocity.

The sharp increase of the wavelength at short times before fingers are visible, as shown in Fig. 14 for $R_{A}=1$ and $\Delta R_{CB}=-1$, might seem surprising but can be explained as follows. Because perturbations are initially dampened, the largest wavelength corresponding to the width $L$ of the system has the largest Fourier amplitude, which increases its weight in the computation of the power-averaged wavelength. This does not happen for larger $\Delta R_{CB}$ as the time when we compute the first $\bar{\lambda}(t)$ is larger than the onset time of the instability.

We now compare the wavelength $\bar{\lambda}$ in the presence and in the absence of reactions. During the linear growth of the instability before merging, $\bar{\lambda}$ is always smaller in the reactive case than for NR case, as already highlighted in Section 3 for specific cases. However, as soon as merging starts, $\bar{\lambda}$ can become larger than its non-reactive counterpart depending on the value of $\Delta R_{CB}$, which affects the merging rate as explained here above.

In conclusion, we have highlighted that the horizontal dynamics characterized by the mean wavelength of the fingering pattern are correlated to the vertical dynamics of the finger tips in the solution. In addition, although the dynamics of the pattern remain similar in all cases, changing the value of $\Delta R_{CB}$ can modify the wavelength in the linear regime and the merging rate in the successive intense and less intense merging regimes. Our results also show that the conclusions drawn from previous studies, i.e. reactions decrease the number of fingers in the pattern 15, 14, are not always true in the fully developed non-linear regime when fingers significantly interact with each other.

4.3 Reaction zone dynamics

Now that we have analyzed the effect of reaction on the convective fingering dynamics, we turn to the influence of convection on the dynamics of reaction fronts. We aim to characterize how and where in the solution the dissolving species A is converted into product C. Recently, we have investigated such reaction dynamics in the absence of convection. To describe the evolution of these reaction-diffusion (RD) fronts, we have derived analytical expressions for the concentration profiles, valid when reaction becomes limited by diffusion, i.e. for times sufficiently larger than the chemical time scale 13, 14, 12. In this limit, the position $z_{f}$ of the reaction front evolves in time as

[TABLE]

where $\eta_{f}$ is a constant depending on the control parameters of the problem. Here, since all species diffuse at the same rate, $\eta_{f}$ depends only on the ratio $\beta$ between the initial concentration $B_{0}$ and the solubility $A_{0}$ as $\eta_{f}=$ erfinv$(1/(1+\beta))$, which is $\approx 0.48$ when $\beta=1$. The RD front delimits two zones: a “reacted” zone above the front, rich in dissolving species A and product C, and an “unreacted” zone below the front, with mostly reactant B as well as some C diffusing towards the bulk solution. In time, the front moves from the interface towards the bulk of the solution as reactant B is progressively depleted. We now examine whether the reaction front $z_{f}$ still progresses as $\sqrt{t}$ in the solution when convection affects the transport dynamics, and evaluate whether modifying $\Delta R_{CB}$ can alter the progression of this reaction front or the evolution of its width.

To quantify the dynamics of the reaction zone in the fully developed non-linear regime, we define the position $z_{f}$ of the reaction front as the first moment of the horizontally-averaged reaction rate profile $\bar{r}(z)$ (see Eq. (21)):

[TABLE]

and the width $w_{f}(t)$ of the reaction front as the width of $\bar{r}(z,t)$ at $\xi=0.1$ of its maximum value, proportional to the second moment of $\bar{r}(z)$ around $z_{f}$:

[TABLE]

After some time, the reaction front position deviates from the diffusive prediction 26, as illustrated in Fig. 15a. For $\Delta R_{CB}\leq 0.2$, including R1 case ($R_{A}=1$, $\Delta R_{CB}=-1$, Figs. 3-5), the reaction front starts to evolve as $t$ and thus moves faster than the RD prediction because of convection, which is coherent with what was observed in horizontal setups where gravity currents occur 45, 46. For $\Delta R_{CB}>0.5$, including R2 case ($R_{A}=1$, $\Delta R_{CB}=1$, Figs. 6-7), the reaction front moves backwards and then stays close to the interface, because amplified convection brings fresh reactant to the interface efficiently. This result is coherent with the concentrations fields in Fig. 7 illustrating that the dissolving species A is consumed as soon as it enters the host solution. Fig. 15b shows that the width $w_{f}$ of the reaction zone increases in time and that this increase varies non monotonically with $\Delta R_{CB}$. The largest values of $w_{f}$ are obtained for intermediate values of $\Delta R_{CB}$ (0 – 0.5). The thinnest reaction zones are thus observed when $\Delta R_{CB}$ is small ($<0$) and the reaction zone cannot extend due to the minimum of density, or when $\Delta R_{CB}$ is large ($>0.5$) and reaction is particularly efficient due to enhanced convective transport.

In summary, we have shown that convection does not always accelerate the progression of the reaction front in the host solution, depending on the value of $\Delta R_{CB}$. When $\Delta R_{CB}$ is large, the reaction takes place in a thin stationary zone close to the interface. When $\Delta R_{CB}$ is small, the reaction zone is also narrow but moves progressively to the bulk at a faster rate than in the absence of convection. The reaction zone is larger for intermediate values of $\Delta R_{CB}$. This means that modifying the composition of the solution, thus impacting $\Delta R_{CB}$, qualitatively affects the reaction zone dynamics during dissolution-driven convection.

5 Storage rate in the presence of reaction

In many applications such as for example CO2 sequestration 4, 5, it is desirable to accelerate the mixing between the dissolving phase and the host phase. Convection increases the mixing between both phases as it increases the movement of the dissolving species A further away from the interface, and increases the flux of A towards the host phase. Moreover, the reaction is expected to also increase the intake of A through consumption. In this context, what are the contributions of convection on the one hand, and of reaction on the other hand, to the evolution of the quantity of dissolved A over time? Is it possible to affect this evolution by selecting given reactants? Does convection coupled with reaction increase the degree of mixing compared with the non-reactive or diffusive-only cases? To answer those questions, we examine the coupled impact of convection and reaction on the evolution of the quantity of dissolved A and more globally on the storage rate of A into the host solution.

5.1 Volume-averaged concentrations

To quantify the storage rate during convective dissolution, we compute the volume-averaged concentration $\langle c_{i}\rangle$ of species $i$ as a function of time as

[TABLE]

From this definition 29, we see that, initially, $\langle A\rangle$ and $\langle C\rangle$ are equal to zero as the dissolving species A and the product C are not present into the host solution yet, while $\langle B\rangle$ is equal to $\beta$. When the host phase is saturated in A, $\langle A\rangle$ is equal to 1 which corresponds to the dimensionless solubility of A in the host solution and thus to its maximum possible concentration. When the reaction is complete, all reactant B has been converted to C as A keeps dissolving into the host phase, so that $\langle C\rangle=\beta$ and $\langle B\rangle=0$.

$\langle A\rangle$ increases in time as species A progressively dissolves into the solution (Fig. 16a). This increase is smaller in the reactive case compared to its non-reactive counterpart, because A is consumed by the reaction. When $\Delta R_{CB}$ is amplified, the increase of $\langle A\rangle$ over time becomes slower; actually, for $\Delta R_{CB}\geq 0.5$, $\langle A\rangle$ is nearly constant for a time period of at least $\approx$ 10000. We explain this steady-state regime in Section 5.2. In addition, as C is produced by the reaction, $\langle C\rangle$ increases and this increase is faster when $\Delta R_{CB}$ increases (Fig. 16b). We define the amount of stored A as $\langle A+C\rangle$, which reflects that species A can be stored in the form of dissolved A or product C. $\langle A+C\rangle$ increases more slowly in the non-reactive case than in the reactive case (Fig. 16c), which means that chemical reactions improve the efficiency of the phase transfer by accelerating the storage process. Most of the stored A is in the form of product C and there is only few dissolved A, as concluded from the comparison of Fig. 16a-c. Further, the larger $\Delta R_{CB}$, the larger is the quantity of A stored as product C.

To better understand why $\langle c_{i}\rangle$ evolve as shown in Fig. 16, we derive evolution equations for these averaged values by integrating over the whole spatial domain the equations (3)-(5) of evolution for the solute concentrations and taking into account the boundary conditions (6)-(9). In the non-reactive case, the quantity of dissolved A in the host phase evolves in time as 29:

[TABLE]

which expresses that $\langle A\rangle$ increases due to the dissolution flux $J$ scaled by the depth of the system $H$, as a deeper system takes more time to achieve saturation. In the presence of reactions, the $\langle c_{i}\rangle$ evolve as

[TABLE]

where $r$ is the volume-averaged reaction rate $\langle AB\rangle$. Equation (31) express that the evolution of the quantities of the solutes depends on the dissolution flux $J$ and the global reaction rate $r$, which we analyze here below.

5.2 Dissolution flux and volume-averaged reaction rate

A dissolution flux through an interface can consist of two different contributions: the convective flux $-A\bm{u}\cdot\bm{n}$, where $\bm{n}$ is the unit vector perpendicular to the interface and oriented towards the outside of the host phase; and the diffusive flux $-\frac{\partial A}{\partial z}$ across the interface. The convective flux across the interface is zero by definition in the theoretical modeling of a monophasic system (see boundary condition 6) and negligeable in biphasic systems 47. We thus compute the interface-averaged dissolution flux $J$ of species A as

[TABLE]

We evaluate the diffusive dissolution flux $J_{D}$ in the non-reactive case as 12

[TABLE]

and the reaction-diffusion flux $J_{RD}$ as 12

[TABLE]

here equal to $2/\sqrt{\pi t}$ as $\beta=1$. Equations (35) and (36) show that for all time $t$, there is a constant ratio $J_{RD}/J_{D}=1+\beta$ between the reactive $J_{RD}$ and non-reactive $J_{D}$ fluxes. Even without convection, the reactive flux is always larger than its non-reactive counterpart as $\beta>0$, and this difference becomes larger when $\beta$ is amplified. Chemical reactions amplify the flux of A across the interface because the consumption of A by the reaction increases the concentration gradient at the origin of the diffusive flux.

As shown in Fig. 17a, the flux initially decreases in time following Eq. (35) (dotted curve) or Eq. (36) (dashed curve) as long as diffusion remains the dominant transport process. After some time, the flux starts to increase because of convection and eventually fluctuates around a steady-state value $J^{*}$. The temporal evolution of $J$ is qualitatively the same for non-reactive and reactive cases. We compute the steady-state flux $J^{*}$ as the average over the last time interval $>3000$ when the variation of the flux with time (least-squares fitted slope) is no more than a small threshold, here arbitrarily chosen as $10^{-6}$. Some fluctuations occur around $J^{*}$ due to reinitiation: the flux increases when protoplumes form as the boundary layer then becomes thinner, and conversely the flux decreases when protoplumes merge with older fingers 9. We compute $J^{*}$ in the non-reactive case as 0.019, in agreement with previous studies 9, 6, 7.

The evolution of the volume-averaged reaction rate $r=\langle AB\rangle$, shown in Fig. 17b, is similar to that of the dissolution flux $J$ (Fig. 17a): when $J$ increases, $r$ increases too. $J$ and $r$ thus appear to be correlated as a more efficient reaction increases the concentration gradient at the origin of the dissolution flux. After a while $r$ also fluctuates around a steady-state value $r^{*}$, which we compute similarly to $J^{*}$ but with a smaller threshold ($4\times 10^{-10}$) for the slope of $r(t)$ given that $r$ is three orders of magnitude smaller than $J$.

We now analyze the variation of the steady-state reaction rate $r^{*}$ (full curve) and scaled flux $J^{*}/H$ (dot-dashed curve) as a function of $\Delta R_{CB}$ in Fig. 18. Both $J^{*}$ and $r^{*}$ increase with $\Delta R_{CB}$, with a change of slope at $\Delta R_{CB}=0$. When $R_{A}=1$, this increase can be described by the empirical fits between $\Delta R_{CB}=-1$ and $1.2$ (see lines in Fig. 18):

[TABLE]

and similarly when $R_{A}=-1$ between $\Delta R_{CB}=0.5$ and $1.2$ (see lines in Fig. 18):

[TABLE]

The increase of $J^{*}$ and $r^{*}$ with $\Delta R_{CB}$, described by Fig. LABEL:eq:_Jr_ra1,eq:_Jr_ra-1, can be explained as follows. When the contribution of the product C to the density increases, convection starts earlier and reaches a larger amplitude (Fig. 13). This increased convection accelerates the transport of fresh reactant B to the interface, which increases the efficiency of the reaction and thus the dissolution flux.

For any value of $\Delta R_{CB}$, the steady-state reaction rate $r^{*}$ is smaller or equal to the scaled steady-state dissolution flux $J^{*}/H$: the reaction is limited by the dissolution rate of A into the host solution. When $R_{A}=1$, we note that if $\Delta R_{CB}<0.5$, $r^{*}$ is typically smaller than $J^{*}/H$, so that A accumulates in the host fluid, and $\langle A\rangle$ increases in time (Eq. (31), Fig. 16a). However, as $r^{*}$ increases faster than $J^{*}/H$ when $\Delta R_{CB}$ increases (see Eqs. (39)-(40)), we see that above $\Delta R_{CB}\geq 0.5$, $r^{*}\approx J^{*}/H$ so that $\langle A\rangle$ remains constant in the steady-state flux regime (Eq. (31), Fig. 16a): as soon as A enters the host fluid, it is consumed by the reaction with B. When $R_{A}=-1$, $J^{*}/H\approx r^{*}$ for all cases as shown in Fig. 18 and expressed by Eqs. (41)-(42).

We can explain the change of slope appearing around $\Delta R_{CB}=0$ in $r^{*}(\Delta R_{CB})$ and $J^{*}(\Delta R_{CB})$ for $R_{A}=1$ (see Eqs. (37)-(40)) as follows. Let us first recall that when all species diffuse at the same rate, the density at the interface is given by $R_{A}+\Delta R_{CB}\beta$, the density at the reaction front is $\Delta R_{CB}\beta$ and the density of the bulk solution is 0 (see Ref. 14 for more details). For $\Delta R_{CB}<0$, the RD density profile in the host solution has a minimum of density. The difference of density at the origin of the instability is then equal to $R_{A}$ and corresponds to the one between the interface and the minimum at the reaction front. Increasing $\Delta R_{CB}$ does not modify this difference but decreases the amplitude of the stabilizing barrier, i.e. the difference of density $-\Delta R_{CB}$ between the minimum and the bulk solution. By contrast, for $\Delta R_{CB}\geq 0$ the density profile in the solution is monotonic. The density difference at the origin of the instability is then between the interface and the bulk solution, and is equal to $\Delta R_{CB}\beta$. This explains why increasing $\Delta R_{CB}$ affects $J^{*}/H$ and $r^{*}$ less when $\Delta R_{CB}<0$ and more when $\Delta R_{CB}\geq 0$.

To analyze the effects of reaction on the steady-state scaled flux $J^{*}/H$, we compare $J^{*}/H$ to its non-reactive counterpart $J^{*}_{NR}/H$ for $R_{A}=1$, plotted in Fig. 18 as a dotted line. For all $\Delta R_{CB}$, $J^{*}/H$ is larger in reactive cases than in NR case. This is not surprising as already without convection, the reaction-diffusion (RD) flux (Eq. (36)) is always $(1+\beta)$ times larger than the diffusive flux (Eq. (35))12. To understand whether this increased steady-state flux is due to RD effects only or also due to convection amplified by reaction, we compare $J^{*}/H$ to the theoretical value of $(1+\beta)J^{*}_{NR}/H$ for the reactive case where only RD effects would affect the dissolution flux. As $J^{*}/H$ increases with $\Delta R_{CB}$, it becomes larger than $(1+\beta)J^{*}_{NR}/H$ for a critical value of $\Delta R_{CB}=0.2$, which corresponds to the critical value computed from the onset time $t_{0}$ (see Fig. 12). In other words, below that critical value, the increase of $J^{*}$ in the reactive case is due to RD effects only, while above that critical value, the increase of $J^{*}$ is also due to convection amplified by reaction.

In summary, the presence of a reactant B in the solution always accelerates the storage of A, even before dissolution-driven convection develops. This can completely change any long-term predictions concerning the fate of A into the host phase. However, depending on the composition of the host solution (impacting $\Delta R_{CB}$), the convection can be amplified by the reaction so that the storage occurs even faster than predicted on the basis of reaction and diffusion alone.

6 Conclusion

We have numerically characterized how the interplay between an A+B$\rightarrow$C reaction and dissolution-driven convection affects the spatio-temporal non-linear dynamics of fingering and reactions fronts as well as the storage efficiency of a phase A into a host fluid phase. Compared to the non-reactive case, reactions can accelerate or slow down fingering development depending on $\Delta R_{CB}$ quantifying the difference between the solutal contributions of product C and of reactant B to density. We have revisited the classification previously established on the basis of a linear stability analysis 13, 14, showing that in some reactive cases convection starts earlier but fingers progress more slowly in the solution than in the non-reactive case, due to non-linear interactions between fingers at later times.

In addition, we have presented different types of convective dynamics, depending on the type of RD density profile building up in the host phase before convection sets in. We have highlighted that when the product C contributes less to the density than the reactant B, a local minimum of density constrains fingering. Hence the reaction front follows the finger tips, but moves downwards at a faster pace than without convection. The fingering pattern originates then mainly from the dissolving species A, progressively accumulating in the host solution as its dissolution takes place at a larger pace than its consumption by the reaction. By contrast, when C contributes much more to the density than B, the reaction, which has a destabilizing effect, takes place mostly in a stationary zone close to the interface. In the steady-state regime, the dissolving species A is consumed as soon as it enters the host solution such that the quantity of A in the host phase remains constant and the fingering pattern is formed mostly by the denser product C. We have also discussed the case where the reaction is at the origin of dissolution-driven convection, because a local maximum forms in the density profile. Due to this specific density profile, fingers then have a characteristic shape and form at a given distance below the interface with phase A.

For all reactive cases, the steady-state flux is larger than its non-reactive counterpart because the consumption of the dissolving species A by the reaction amplifies the concentration gradient at the origin of the dissolution flux. Both steady-state flux and reaction rate increase with $\Delta R_{CB}$ as convection develops earlier and becomes more intense, amplifying the mixing between the phase A and the host phase. Although this model could be extended to include differential diffusivity effects 12 and a variable solubility depending on the solute concentrations 42, our results already highlight that selecting an appropriate composition for the host phase allows to maximise the positive effect of reaction for amplifying the storage rate of A into the host fluid.

These conclusions are useful to predict the fate of CO2 during its sequestration in subsurface formations and to select storage sites with geochemical reactions optimal in enhancing convective dissolution. Knowing the kinetic properties of the geochemical reaction is not enough to quantify the storage rate; one also needs to know the contributions to density of all dissolved species. For other applications where convection enhances mass transfer, controlling the properties of the dissolution-driven convection should become possible by selecting the appropriate reactant to be dissolved in the host solution.

Acknowledgments

Funding by PRODEX, ARC CONVINCE, ARC PIONEER, MIS-FNRS PYRAMID and PDR-FNRS FORECAST projects is gratefully acknowledged. We thank V. Moureau and G. Lartigue for their training and support concerning YALES2 code, as well as F. Brau, C. Rana and V. Upadhyay for useful comments and scientific discussions on this study.