Nagatsu and De Wit [“Viscous fingering of a miscible reactive A + B → C interface for an infinitely fast chemical reaction: Nonlinear simulations,” Phys. Fluids **23**, 043103 (2011)] simulated the nonlinear evolution of reactive miscible viscous fingering (VF) where the viscosity of the more viscous displaced fluid was changed by an instantaneous A + B → C chemical reaction. They analyzed the dynamics from the viewpoint of an underlying viscosity profile reconstructed from the concentrations of chemical species obtained by a one-dimensional diffusion–reaction equation. The present study develops a mathematical model for reactive miscible VF where the viscosity of the less viscous displacing fluid is changed by an instantaneous A + B → C chemical reaction. We obtain the same underlying viscosity profile as Nagatsu and De Wit by employing appropriate parameters. We perform numerical simulations of the nonlinear evolution of VF under these appropriate parameters. The results show that the present numerical solutions are exactly the same as those obtained by Nagatsu and De Wit, i.e., the same VF pattern is obtained. This numerically proves that the effects of a viscosity change by the chemical reaction on VF are independent of whether the viscosity of the displaced or displacing fluid changes. We obtain a mathematical formula to describe the switch from the parameters used by Nagatsu and De Wit to those used in this study to obtain the same shape of the underlying viscosity profile. This finding will lead to easier and more flexible VF chemical control methods in geoscience processes, where it is difficult to manipulate the properties of more viscous fluids, by manipulating the properties of less viscous fluids.

## I. INTRODUCTION

When a more viscous fluid is displaced by a less viscous one in porous media or a Hele-Shaw cell, the interface between the two fluids becomes hydrodynamically unstable and a finger-like pattern is formed. This phenomenon is called viscous fingering (VF).^{1–3} VF instability can occur in a variety of processes, including enhanced oil recovery (EOR),^{4} CO_{2} sequestration,^{5} and chromatographic separation processes.^{6} Therefore, various studies have attempted to elucidate the basic characteristics of VF.^{7,8} Classically, VF is divided into two categories, namely, immiscible VF and miscible VF. The control parameters governing the dynamics of immiscible VF and miscible VF are the capillary number Ca and the Péclet number Pe, respectively. Here, Ca is the ratio between the viscous and interfacial forces, whereas Pe is the ratio between the mass transfer rate by convection and that by diffusion. Note that the viscosity ratio of the more and less viscous fluids primarily governs the dynamics. For both systems, nonlinear propagation of VF is governed by different mechanisms of shielding, spreading, and splitting. Shielding is the phenomenon whereby a finger that is slightly ahead of its neighboring fingers quickly outruns them and shields them from further growth. Spreading and splitting are the phenomena in which a finger that spreads until it reaches a certain width becomes unstable and splits.^{7}

In processes such as EOR,^{9} CO_{2} sequestration,^{10} and medical applications,^{11} VF is often accompanied by a chemical reaction. Experimental and theoretical studies on coupling between VF and chemical reactions have recently been conducted,^{12–41} and several review articles on this issue are available.^{42–44} Here, we are especially interested in miscible VF with an A + B → C chemical reaction that changes the viscosity.^{17,20,22–28,30–32,34,39,41} Nagatsu *et al.*^{17} conducted the first experiments on reactive miscible VF in a radial Hele-Shaw cell and found that an A + B → C type of reaction actively changes the viscosity and, hence, affects VF. In these experiments, more viscous solutions of polymers (the viscosity of which varied with pH) were displaced by less viscous miscible solutions of acid or base. The reaction was very fast. They found that an increase in viscosity by the reaction induced wider fingers and suppressed the shielding effect, producing a pattern covering a larger area of the displacing solution. Conversely, a decrease in viscosity led to thinner fingers with a stronger shielding effect, producing a pattern covering a smaller area of the displacing solution. A mathematical model has been developed to describe the experimental results of Nagatsu *et al.*^{23,26,28,41} These theoretical studies considered a solution including reactant A that pushed another solution including reactant B with an A + B → C chemical reaction, where the solution’s viscosity depended on the concentrations of B and C. In addition, they have attempted to classify the dynamics through the underlying viscosity profiles reconstructed from the concentration profiles of chemical species obtained from one-dimensional diffusion–reaction equations. In these studies, Nagatsu and De Wit^{28} performed numerical simulations of the VF nonlinear dynamics under the condition of infinite D_{a} in rectilinear geometry (D_{a} is defined as the ratio between the characteristic time of fluid motion and that of the chemical reaction; infinite D_{a} denotes a reaction rate that is instantaneously fast). They numerically showed that the shielding effect was suppressed and a denser fingering pattern was obtained when the viscosity was increased, whereas the shielding effect was enhanced and a less-dense fingering pattern was obtained when the viscosity was decreased. The results were in good agreement with the experimental results of Nagatsu *et al.*^{17}

Nagatsu *et al.*^{24} performed a VF experiment in which the viscosity of a less viscous fluid was changed by chemical reactions. They used a polymer solution whose viscosity depends on pH as the less viscous fluid and glycerol or starch syrup solution including reactants (acid or base) as the more viscous fluid. They reported similar results to those obtained by Nagatsu *et al.*^{17} Namely, finger widening and suppression of the shielding effects occurred when the viscosity was increased. On the contrary, finger thinning and enhancement of the shielding effect occurred when the viscosity was decreased. These experimental results show that the effects of the viscosity changes on the VF pattern are equivalent, regardless of whether the viscosity of the less or more viscous fluid is changed. They explained that the one-dimensional underlying viscosity profile became equivalent for both cases, and this was the cause of the equivalent influences of viscosity changes of the more and less viscous fluids on the nonlinear evolution of VF. However, to complete the investigation of this problem, numerical simulations are necessary. In this context, the present study develops a mathematical model for reactive miscible VF in which the viscosity of the less viscous displacing fluid is changed by an A + B → C chemical reaction. Here, we consider that a less viscous solution containing reactant A displaces a more viscous solution containing reactant B and a highly viscous solvent D, in which the viscosity depends on the concentrations of A, C, and D. We investigate the one-dimensional underlying viscosity profile described above and perform numerical simulations of the nonlinear evolution of reactive miscible VF. We emphasize that the viscosity depended on the concentrations of chemical species B and C; in other words, the viscosity of the displaced more viscous fluid was changed, in other studies except that of Nagatsu *et al.*,^{24} which experimentally and theoretically investigated miscible VF with the A + B → C reaction changing the viscosity.^{17,20,22,23,25–28,31,32,34,39,41}

## II. MODEL

Figure 1 shows the model system and the initial concentration conditions. We consider a uniform porous medium of length $Lx\u2032$, $Ly\u2032$ and permeability *κ*. A less viscous fluid with a viscosity *μ* = *μ*_{A} containing reactant A pushes a more viscous fluid that contains reactant B and a highly viscous solvent D with a viscosity *μ*_{D}, from left to right with the flow velocity *U*. Reactants A and B are present at the same initial concentrations, *a*_{0} = *b*_{0}, respectively. Once reactants A and B contact at the interface, a simple A + B → C chemical reaction occurs immediately, generating product C, which has a viscosity of *μ*_{C}. The highly viscous solvent D, with an initial concentration of *d*_{0}, is inert to this chemical reaction (glycerol and starch syrup used by Nagatsu *et al.*^{24}). As in the study of Nagatsu and De Wit,^{28} the present study added a chemical species *E* to the less viscous fluid, corresponding to the dye and represented by the initial concentration *e*_{0}.

This system is considered to be incompressible and neutrally buoyant. The governing equations are the continuity equation (1), momentum conservation equation (2), and chemical species conservation equations (3)–(7), which can be written, respectively, as

Here, ** u** = (

*u*,

*v*) represents the two-dimensional flow velocity and

*p*is the pressure. The parameters

*a*,

*b*,

*c*,

*d*, and

*e*represent the concentrations of chemical species A, B, C, D, and E, respectively, which are molar concentrations and are expressed in dimensions of moles/

*L*

^{3}(

*L*means the dimension of length). As the mass concentration of chemical species is frequently employed in chemically reactive fluid dynamics, the derivation of Eqs. (3)–(7) is given in the Appendix.

*k*is a reaction rate constant, and

*D*

_{A,B,C,D,E}represents the diffusion coefficient of A, B, C, D, and E, respectively. Here, we assume

*D*

_{A}=

*D*

_{B}=

*D*

_{C}=

*D*

_{D}=

*D*

_{E}=

*D*.

We define the viscosities of solutions in which species A, B, C, and D are dissolved at a concentration of *f*_{0} as *μ*_{A}, *μ*_{B}, *μ*_{C}, and *μ*_{D}. The viscosity of the solution is given by the following expression, where a component F with viscosity *μ* = *μ*_{0} is used as a reference:

Here, *R*_{a}, *R*_{c}, *R*_{d} are natural logarithmic ratios of the viscosity defined by

This definition is different from that employed in the previous studies, where the viscosity depends on the concentrations of B and C and the viscosities of solutions in which species A, B, and C are dissolved at a concentration of *a*_{0} are *μ*_{A}, *μ*_{B}, and *μ*_{C}.^{23,25,26,28,30,31,41} As the viscosity depends on *a* in the present study, it is necessary to normalize the viscosity with another component whose viscosity does not depend on the concentration.

First, we switch to a reference frame moving with the injection speed *U*.^{45} We nondimensionalize the obtained equations using the characteristic velocity *U*, the hydrodynamical time *τ*_{h} = *D*/*U*^{2}, and the length *L*_{h} = *D*/*U*. Viscosity is normalized by *μ*_{0} and pressure by *μ*_{0}*D*/κ. The concentrations of species A, B, C, and D are scaled by *f*_{0} and the concentration of dye E by *e*_{0}. The width and length of the system, $Lx=ULx\u2032/D$ and $Ly=ULy\u2032/D$, respectively, are dimensionless.^{45} Taking the curl of Darcy’s law, and introducing the stream function *ψ*(*x*, *y*) defined by *u* = *∂ψ*/*∂y* and *v* = −*∂ψ*/*∂x*, yields the following set of dimensionless equations:

Here, for simplicity, no symbols are used to represent dimensionless numbers. The subscripts *x*, *y*, and *t* represent the derivatives with respect to *x*, *y*, and *t*, respectively. *D*_{a} = *Dkf*_{0}/*U*^{2} is a dimensionless number given by the ratio of the reaction rate to the flow rate, known as the Damköhler number.

Following the procedure described by Nagatsu and De Wit,^{28} Eqs. (10)–(15) are solved under the condition *D*_{a} → ∞, that is, the condition of an instantaneous reaction. Subtracting Eq. (11) from Eq. (12) and adding Eq. (11) to Eq. (13), we get, respectively,

We then introduce two conserved scalars *Z*_{1} = *a* − *b* and *Z*_{2} = *a* + *c*, which are independent of the chemical reaction. These conserved scalars can be normalized in the form

where *i* = 1, 2, while *Z*_{i,A0} and *Z*_{i,B0} are the values of *Z*_{i} in the pure displacing and displaced fluid, respectively. Therefore in the displacing fluid, *z*_{i} = 1, while in the displaced fluid, *z*_{i} = 0. Regarding *Z*_{1} = *a* − *b*, because *Z*_{1,A0} = 1 and *Z*_{1,B0} = −1, we obtain

Regarding *Z*_{2} = *a* + *c*, because *Z*_{2,A0} = 1 and *Z*_{2,B0} = 0, we obtain

It should be noted that now *Z*_{1} and *Z*_{2} are normalized between 0 and 1 and thus *z*_{1} = *z*_{2} = *z*. Under the condition of infinite *D*_{a}, the chemical reaction takes place in an infinitesimally thin region called here the reaction front. At the reaction front, the concentration of reactants *a* and *b* is 0, while *z* becomes *z*_{st} (the subscript st stands for stoichiometry) that takes the value *z*_{st} = 0.5 from Eq. (20). When *z* ≤ *z*_{st}, *a* = 0, *b* = 1 − 2*z* from Eq. (20) and *c* = *z* from Eq. (21). When *z* ≥ *z*_{st}, *b* = 0, *a* = −1 + 2*z* from Eq. (20) and *c* = 1 − *z* from Eq. (21). In summary, under the condition *D*_{a} → ∞, the governing equations to be solved can be reduced to the following:

along with the condition that for *z* ≤ *z*_{st}, we have

whereas for *z* ≥ *z*_{st},

Note that the dye concentration *e* follows exactly the same dynamics as *z*.

In solving Eqs. (22)–(26), the pseudospectral method was used with periodic boundary conditions in the *x* and *y* directions. The size of the system was *L*_{x} × *L*_{y} = 8192 × 2048. The initial conditions of the stream function and the product concentration were *ψ*(*x*, *y*) = 0 and *c*(*x*, *y*) = 0 for all (*x*, *y*), respectively. For the concentrations, we considered a step front between the reactive solutions of A and B with *c* = 0 everywhere and added random noise in the front to trigger the instability within a reasonable computing time. Here, the initial concentration of each chemical species was fixed, that is, *a*_{0} = *b*_{0} = *d*_{0} = *e*_{0} = 1. In this study, we consider a viscously unstable system in the absence of the reaction, which requires *R*_{a} < *R*_{d}. Here, we set *R*_{a} = 1 and *R*_{d} = 3, −3 ≤ *R*_{c} ≤ 5, to investigate the effect of the viscosity of species C generated by the chemical reaction.

## III. UNDERLYING VISCOSITY PROFILE

Before analyzing the results of the nonlinear simulations, we first investigate the underlying viscosity profile. The concentration profiles of the chemical species given by one-dimensional reaction–diffusion equations for both reactive and nonreactive cases [which are obtained by setting *ψ*_{x} = *ψ*_{y} = 0 in Eqs. (10)–(14)] and the viscosity profiles reconstructed by formula (16) are shown in Fig. 2. The horizontal axis is the self-similar variable $\eta =x/4t$. In the absence of a reaction, the concentration profiles of the passive scalars A, B, and D evolve following an error function [Fig. 2(a)], which gives an unstable viscosity profile increasing monotonically with *η* because *R*_{a} < *R*_{d}. In the presence of a chemical reaction, the profile of reactants A and B changes and product C appears. However, the profile of D remains unchanged because it is inert to the reaction. In the special case where *R*_{c} = 1 (=*R*_{a}), the viscosity profile equals the nonreactive one because the consumption of A by the reaction is exactly balanced in the viscosity profile by the production of C, which has the same viscosity as A. When *R*_{c} > 1 (=*R*_{a}), product C is more viscous than A and the viscosity increases in the reaction zone. The viscosity profile becomes nonmonotonic with a maximum once *R*_{c} > 3. Analogous but symmetrical situations are obtained for *R*_{c} < 1 (=*R*_{a}), where nonmonotonic profiles with a minimum develop if *R*_{c} < −1. Note in Fig. 2(b) that the two conditions corresponding to a fixed |*R*_{c} − *R*_{a}| yield symmetric viscosity profiles with respect to the nonreactive case *R*_{c} = *R*_{a}. We emphasize that the shapes of the viscosity profiles shown in Fig. 2(b) are exactly the same as those in Fig. 5(b) in the study of Nagatsu and De Wit,^{28} although the range of ln(*μ*) is different [0 ≤ ln(*μ*) ≤ 4 in the present study, whereas −1 ≤ ln(*μ*) ≤ 3 in the study of Nagatsu and De Wit^{28}].

## IV. RESULTS AND DISCUSSION

Figures 3 and 4 show the results of the numerical simulations for viscosity increase cases (*R*_{c} ≥ 1) and viscosity decrease cases (*R*_{c} ≤ 1), respectively. These simulations were performed under the same initial disturbance, time, and spatial steps as those employed by Nagatsu and De Wit.^{28} Here, we show the temporal evolution of the concentration field of dye E. The gray scale ranges from 0 in white to 1 in black. The results in Figs. 3 and 4 are exactly the same as those in Figs. 6 and 7 in the study of Nagatsu and De Wit.^{28} This shows that exactly the same VF pattern for dye is obtained when the underlying viscosity profile is equivalent. This also numerically demonstrates that the effects of viscosity change by the reaction are identical, regardless of whether the viscosity of the less viscous displacing fluid or more viscous displaced fluid is changed, as experimentally shown by Nagatsu *et al.*^{24} From these results, we can conclude that the independence of the viscosity change effects on reactive miscible VF on which fluid causes the viscosity change arises from the fact that the underlying viscosity profile in the reactive zone changes in a similar way.

As mentioned above, because the numerical results for the concentration fields of the dye shown in Figs. 3 and 4 are exactly the same as those shown in Figs. 6 and 7 in the study of Nagatsu and De Wit,^{28} the results for the finger length and finger density shown in Figs. 9 and 10 in the study of Nagatsu and De Wit^{28} are also applicable to this study. In other words, when the graph legends from *R*_{c} = −2 to 6 in the study of Nagatsu and De Wit^{28} are replaced by *R*_{c} = −3 to 5, respectively, the same results are obtained in the present study.

Figure 5 shows the concentration fields of A, B, C, D, and dye E and viscosity and velocity fields for various *R*_{c}. This figure corresponds to Fig. 12 in the study of Nagatsu and De Wit.^{28} We emphasize that Fig. 5 includes the concentration field of D. Note that we show the velocity in the *x*-direction, *u*, as the velocity field, although the magnitude of velocity $u2+v2$ was shown in the study of Nagatsu and De Wit.^{28} In Fig. 5, the gray scale ranges from 0 in white to 1 in black for the concentration fields of A, B, D, and E, from 0 in white to 0.5 in black for the concentration field of C, and from −0.7 in white to 2.5 in black for the velocity field. For the velocity fields, the gray scale ranges from the minimum value in white to the maximum value in black in each *R*_{c} condition. The concentration fields of A, B, C, and E and the viscosity fields are exactly the same as those in the study of Nagatsu and De Wit.^{28} This is considered to be because the concentration profiles of A, B, and C [shown in Fig. 2(a) in the present study] are exactly the same as those in Fig. 5(a) in the study of Nagatsu and De Wit.^{28} The complete equivalence of dye and viscosity fields between the present study and the study of Nagatsu and De Wit^{28} arises from the complete equivalence of the underlying one-dimensional viscosity profiles. We have found that the concentration fields of B and D are different, with the concentration of B being less than that of D in the fingers, although they contain the same initial non-dimensional concentrations. This is clearly caused by the chemical reaction in which B is consumed in the finger. We note that as mentioned by Nagatsu and De Wit,^{28} the concentration of A is lower than that of E in the finger, although they have the same initial dimensionless concentrations; this is caused by the chemical reaction. Additionally, the velocity fields are essentially the same as those in the study of Nagatsu and De Wit.^{28} We have confirmed that the velocity field shown in Fig. 5 indicates that larger-velocity regions develop in adjacent fingers in the dye representation for the viscosity increase case [(b) and (c)], which suppresses the shielding effect, whereas these regions develop in the most advanced finger for the viscosity decrease case [(d) and (e)], which enhances the shielding effect. The detailed mechanism for the suppression or enhancement of the shielding effect was described by Nagatsu and De Wit.^{28}

## V. MATHEMATICAL STUDY

As mentioned above, the shape of the underlying one-dimensional viscosity profile is equivalent, regardless of whether the viscosity of the less or more viscous fluid is changed by the reaction. Note that that the complete coincidence of the underlying one-dimensional viscosity profile is not necessary. In fact, the range of ln(*μ*) profiles in this study is different from that obtained by Nagatsu and De Wit.^{28} In this section, we derive a mathematical formula that explains the switch from the parameters used by Nagatsu and De Wit^{28} to those used in this study, to achieve the same shape of the underlying one-dimensional viscosity profile.

In this study and that of Nagatsu and De Wit,^{28} analytical solutions are realized for one-dimensional concentrations of chemical species by specifying an infinite value of D_{a}. The following equations provide the solutions for the non-reactive and reactive cases, respectively, which are illustrated in Fig. 2:

Nagatsu and De Wit^{28} did not consider the solutions for the *d* species. They represented the underlying viscosity profile using the following equation:

By substituting Eqs. (27) and (28) into Eq. (29), we obtain the analytical solutions of the one-dimensional underlying viscosity profile, wherein the viscosity depends on *b* and *c*, shown in Fig. 5 in the study of Nagatsu and De Wit,^{28}

Similarly, by substituting Eqs. (27) and (28) into Eq. (16), we obtain the analytical solutions of the one-dimensional underlying viscosity profile, wherein the viscosity depends on *a*, *c*, and *d*, shown in Fig. 2,

The condition required for coincidence of the shape of the viscosity profile given by Eqs. (30) and (31) with that given by Eqs. (32) and (33) is the equivalence of the gradients (derivative with respect to *η*) for the non-reactive case, the reactive case with *η* ≤ 0, and the reactive case with *η* ≥ 0. That is,

where *R*_{c1} is *R*_{c} in the study of Nagatsu and De Wit,^{28} whereas *R*_{c2} is *R*_{c} in this study. From Eq. (34), the following relations are obtained:

In the study of Nagatsu and De Wit,^{28}^{,}*R*_{b} = 2. Then, *R*_{a} is set to 1 to satisfy Eq. (35); *R*_{d} = 3 and *R*_{c2} = *R*_{c1} − 1 are obtained. Actually, under the condition *R*_{c2} = *R*_{c1} − 1, the same shape of the viscosity profile is obtained as in the study of Nagatsu and De Wit.^{28} If we set *R*_{b} = 3 and *R*_{a} = 0.5, *R*_{d} = 3.5 and *R*_{c2} = *R*_{c1} − 2.5 are obtained for Eq. (35). To verify this calculation, the underlying viscosity profiles under the condition of *R*_{b} = 3, where the viscosity depends on *b* and *c*, and of *R*_{a} = 0.5, where the viscosity depends on *a*, *c*, and *d*, are shown in Fig. 6. We confirm that the conditions of *R*_{d} = 3.5 and *R*_{c2} = *R*_{c1} − 2.5 give the same shape of the profile although the range of ln(*μ*) is different. We emphasize that the formulations represented by Eq. (35) for switching from the parameters used in the study of Nagatsu and De Wit^{28} to those used in this study to obtain the same numerical results can be applicable under the condition where the reactant concentrations are identical and diffusion coefficients of all chemical species are the same.

Herein, we discuss how the same underlying viscosity profiles produce the same numerical simulation results. In the study of Nagatsu and De Wit (2011), the vorticity formulation is given as

where we denote the spatial derivative of the error function as *d*_{e} and we use the following equation:

A similar analysis is done for the present system, in which the vorticity formulation is given as Eq. (10). By inserting Eqs. (27) and (28) into Eq. (10), we obtain the following equations for the non-reactive and reactive cases, respectively:

After equating Eq. (37) to Eq. (40) and Eq. (39) to Eq. (41), we obtain Eq. (34). This originates from the vorticity equation; thus, it can explain that satisfaction of Eq. (34) provides the same numerical results of Eqs. (10) and (36) in a more direct manner. In addition, this discussion can explain that the condition of *R*_{b} = *R*_{c} provides the equivalent vorticity equation for all *η* in the non-reactive and reactive cases [Eqs. (37) and (39)].

## VI. CONCLUSIONS

The present study has developed a mathematical model of reactive miscible VF in which the viscosity of the less viscous displacing fluid is changed by an instantaneous A + B → C chemical reaction. Here, a less viscous solution containing reactant A displaces a more viscous solution containing reactant B and a highly viscous solvent D, in which the viscosity depends on the concentrations of A, C, and D. We obtained the same underlying viscosity profile as that in the study of Nagatsu and De Wit^{28} by employing appropriate parameters. The results of numerical simulations of the nonlinear evolution of VF for these parameters show that the present numerical solutions are exactly the same as those obtained by Nagatsu and De Wit,^{28} i.e., the same VF pattern is obtained. This result shows that the same VF pattern for dye is obtained when the underlying viscosity profile is equivalent. Thus, we have numerically demonstrated that the effects of the viscosity change by the reaction are identical regardless of whether the viscosity of the less viscous displacing fluid or more viscous displaced fluid is changed, as experimentally shown by Nagatsu *et al.*,^{24} and that the origin of these effects is that the underlying viscosity profile in the reactive zone changes in a similar way. In addition, we obtained a mathematical formula to describe the switch from the parameters used in the study by Nagatsu and De Wit^{28} to those in this study to obtain the same shape of the underlying viscosity profile, expressed by Eq. (35). Furthermore, we discussed our observations regarding the same underlying viscosity profiles producing the same numerical simulation results.

From the viewpoint of VF control, the independence of the effects of the reaction on the VF pattern on which fluid changes viscosity shows that the same effect can be obtained by operating a less viscous or a more viscous fluid. This is important for effective chemical control of VF in geoscience processes such as EOR and CO_{2} sequestration because it is difficult to manipulate the properties of more viscous fluids in these processes. The results presented in this paper will lead to easier and more flexible VF chemical control methods by manipulating the properties of less viscous fluids.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

We thank Professor A. De Wit of Université Libre de Bruxelles for providing us with the numerical code used in the study of Nagatsu and De Wit.^{28} We would like to thank Editage (www.editage.com) for English language editing.

### APPENDIX: DERIVATION OF EQS. (3)–(7)

Here, the derivation of Eqs. (3)–(5) is shown. First, the mass conservation of chemical species A, B, and C in a diluted liquid solution where the mass concentration of the mixture, *ρ*, which has the dimension of *M*/*L*^{3} (*M* is the dimension of mass), and diffusion coefficients of chemical species A, B, and C, *D*_{A,B,C}, are constant is written as^{46}

where *ρ*_{A}, *ρ*_{B}, *ρ*_{c} are the mass concentration of chemical species A, B, and C, whose dimensions are *M*/*L*^{3}. In addition, *r*_{A} and *r*_{B} are the mass consumption rates of chemical species A and B, which are the reactants, while *r*_{C} is the mass production rate of chemical species C, the product. The dimensions of *r*_{A}, *r*_{B}, and *r*_{C} are *M*/(*L*^{3}*t*), where *t* is the dimension of time. When Eqs. (A1)–(A3) are divided by the molecular weight of chemical species A, B, and C, *M*_{A}, *M*_{B}, and *M*_{C}, whose dimensions are *M*/mole, the following equations are obtained:

Here, *a* = *ρ*_{A}/*M*_{A}, *b* = *ρ*_{B}/*M*_{B}, and *c* = *ρ*_{C}/*M*_{C}. The reaction terms are now expressed with the dimension of moles/(*L*^{3}*t*) as follows:

By combining Eqs. (A4)–(A6) with Eq. (A7), Eqs. (3)–(5), as represented in the main text, are obtained. From Eq. (A7), we obtain

where *M*_{A} + *M*_{B} = *M*_{C} because the equimolar reaction A + B → C is considered. We note that the total mass of chemical species A, B, and C is conserved through the reaction, whereas the total number of moles of chemical species A, B, and C is not. Similarly, Eqs. (6) and (7) in the main text, without the reaction term, are obtained.

## REFERENCES

_{2}in deep saline aquifers: Analytical solution for CO

_{2}plume evolution during injection

_{2}sequestration in a radial Hele–Shaw cell via an interfacial chemical reaction