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.

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 CO2 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 CO2 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 Wit28 performed numerical simulations of the VF nonlinear dynamics under the condition of infinite Da in rectilinear geometry (Da is defined as the ratio between the characteristic time of fluid motion and that of the chemical reaction; infinite Da 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

Figure 1 shows the model system and the initial concentration conditions. We consider a uniform porous medium of length Lx, Ly 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, a0 = b0, 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 d0, 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 e0.

FIG. 1.

Sketch of the system.

FIG. 1.

Sketch of the system.

Close modal

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

u=0,
(1)
p=μ(a,c,d)κu,
(2)
at+ua=DA2akab,
(3)
bt+ub=DB2bkab,
(4)
ct+uc=DC2c+kab,
(5)
dt+ud=DD2d,
(6)
et+ue=DE2e.
(7)

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/L3 (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 DA,B,C,D,E represents the diffusion coefficient of A, B, C, D, and E, respectively. Here, we assume DA = DB = DC = DD = DE = D.

We define the viscosities of solutions in which species A, B, C, and D are dissolved at a concentration of f0 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:

μa,c,d=μ0expRaaf0+Rccf0+Rddf0.
(8)

Here, Ra, Rc, Rd are natural logarithmic ratios of the viscosity defined by

Ra=lnμAμ0,Rc=lnμCμ0,Rd=lnμDμ0.
(9)

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 a0 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/U2, and the length Lh = D/U. Viscosity is normalized by μ0 and pressure by μ0D/κ. The concentrations of species A, B, C, and D are scaled by f0 and the concentration of dye E by e0. The width and length of the system, Lx=ULx/D and Ly=ULy/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:

2ψ=Raψxax+ψyay+ay+Rcψxcx+ψycy+cy+Rdψxdx+ψydy+dy,
(10)
at+axψyayψx=2aDaab,
(11)
bt+bxψybyψx=2bDaab,
(12)
ct+cxψycyψx=2c+Daab,
(13)
dt+dxψydyψx=2d,
(14)
et+exψyeyψx=2e,
(15)
μa,c,d=expRaa+Rcc+Rdd.
(16)

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. Da = Dkf0/U2 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 Da → ∞, that is, the condition of an instantaneous reaction. Subtracting Eq. (11) from Eq. (12) and adding Eq. (11) to Eq. (13), we get, respectively,

abt+abxψyabyψx=2ab,
(17)
a+ct+a+cxψya+cyψx=2a+c.
(18)

We then introduce two conserved scalars Z1 = ab and Z2 = a + c, which are independent of the chemical reaction. These conserved scalars can be normalized in the form

zi=ZiZi,B0Zi,A0Zi,B0,
(19)

where i = 1, 2, while Zi,A0 and Zi,B0 are the values of Zi in the pure displacing and displaced fluid, respectively. Therefore in the displacing fluid, zi = 1, while in the displaced fluid, zi = 0. Regarding Z1 = ab, because Z1,A0 = 1 and Z1,B0 = −1, we obtain

z1=ab+12.
(20)

Regarding Z2 = a + c, because Z2,A0 = 1 and Z2,B0 = 0, we obtain

z2=a+c.
(21)

It should be noted that now Z1 and Z2 are normalized between 0 and 1 and thus z1 = z2 = z. Under the condition of infinite Da, 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 zst (the subscript st stands for stoichiometry) that takes the value zst = 0.5 from Eq. (20). When zzst, a = 0, b = 1 − 2z from Eq. (20) and c = z from Eq. (21). When zzst, b = 0, a = −1 + 2z from Eq. (20) and c = 1 − z from Eq. (21). In summary, under the condition Da → ∞, the governing equations to be solved can be reduced to the following:

2ψ=Raψxax+ψyay+ay+Rcψxcx+ψycy+cy+Rdψxdx+ψydy+dy,
(22)
zt+zxψyzyψx=2z,
(23)
dt+dxψydyψx=2d,
(24)

along with the condition that for zzst, we have

(a,b,c)=(0,12z,z),
(25)

whereas for zzst,

(a,b,c)=(1+2z,0,1z).
(26)

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 Lx × Ly = 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, a0 = b0 = d0 = e0 = 1. In this study, we consider a viscously unstable system in the absence of the reaction, which requires Ra < Rd. Here, we set Ra = 1 and Rd = 3, −3 ≤ Rc ≤ 5, to investigate the effect of the viscosity of species C generated by the chemical reaction.

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 η=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 Ra < Rd. 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 Rc = 1 (=Ra), 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 Rc > 1 (=Ra), product C is more viscous than A and the viscosity increases in the reaction zone. The viscosity profile becomes nonmonotonic with a maximum once Rc > 3. Analogous but symmetrical situations are obtained for Rc < 1 (=Ra), where nonmonotonic profiles with a minimum develop if Rc < −1. Note in Fig. 2(b) that the two conditions corresponding to a fixed |RcRa| yield symmetric viscosity profiles with respect to the nonreactive case Rc = Ra. 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 Wit28].

FIG. 2.

1D reaction–diffusion: (a) dimensionless concentration profiles of the various chemical species and (b) viscosity profiles reconstructed using Eq. (16) for various Rc at Ra = 1, Rd = 3. The notation “NR” stands for “nonreactive.”

FIG. 2.

1D reaction–diffusion: (a) dimensionless concentration profiles of the various chemical species and (b) viscosity profiles reconstructed using Eq. (16) for various Rc at Ra = 1, Rd = 3. The notation “NR” stands for “nonreactive.”

Close modal

Figures 3 and 4 show the results of the numerical simulations for viscosity increase cases (Rc ≥ 1) and viscosity decrease cases (Rc ≤ 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.

FIG. 3.

VF dynamics of dye E for various RcRa = 1 at Ra = 1, Rd = 3 at time t = 400, 800, 1600, 2400, and 3200 (from top to bottom). The red lines represent the position of the initial interface between the less and more viscous fluids.

FIG. 3.

VF dynamics of dye E for various RcRa = 1 at Ra = 1, Rd = 3 at time t = 400, 800, 1600, 2400, and 3200 (from top to bottom). The red lines represent the position of the initial interface between the less and more viscous fluids.

Close modal
FIG. 4.

VF dynamics of dye E for various RcRa = 1 at Ra = 1, Rd = 3 at time t = 400, 800, 1600, and 2400 (from top to bottom). The red lines represent the position of the initial interface between the less and more viscous fluids.

FIG. 4.

VF dynamics of dye E for various RcRa = 1 at Ra = 1, Rd = 3 at time t = 400, 800, 1600, and 2400 (from top to bottom). The red lines represent the position of the initial interface between the less and more viscous fluids.

Close modal

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 Wit28 are also applicable to this study. In other words, when the graph legends from Rc = −2 to 6 in the study of Nagatsu and De Wit28 are replaced by Rc = −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 Rc. 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 Rc 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 Wit28 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 

FIG. 5.

Concentration of A, B, C, D, and dye E and viscosity and velocity fields shown from top to bottom for Ra = 1, Rd = 3, and (a) Rc = 1 at time t = 2000, (b) Rc = 3 at time t = 2000, (c) Rc = 5 at time t = 2800, (d) Rc = −1 at time t = 2000, and (e) Rc = −3 at time t = 2000.

FIG. 5.

Concentration of A, B, C, D, and dye E and viscosity and velocity fields shown from top to bottom for Ra = 1, Rd = 3, and (a) Rc = 1 at time t = 2000, (b) Rc = 3 at time t = 2000, (c) Rc = 5 at time t = 2800, (d) Rc = −1 at time t = 2000, and (e) Rc = −3 at time t = 2000.

Close modal

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 Wit28 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 Da. The following equations provide the solutions for the non-reactive and reactive cases, respectively, which are illustrated in Fig. 2:

a=1erf(η)2,b=1+erf(η)2,c=0,d=1+erf(η)2, for all η,
(27)
a,b,c,d=erfη,0,1+erfη2,1+erf(η)2,η<00,erfη,1erfη2,1+erf(η)2,η0.
(28)

Nagatsu and De Wit28 did not consider the solutions for the d species. They represented the underlying viscosity profile using the following equation:

lnμb,c=Rbb+Rcc.
(29)

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 

lnμ=Rb+Rberf(η)2 for all η,
(30)
lnμ=Rc+Rcerf(η)2,η<0Rc+(2RbRc)erf(η)2,η0.
(31)

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,

μ=Ra+Rd+RdRaerf(η)2 for all η,
(32)
μ=Rc+Rd+Rc+Rd2Raerfη2,η<0Rc+Rd+(RdRc)erf(η)2,η0.
(33)

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,

Rb2=RdRa2,Rc12=Rc2+Rd2Ra2,2RbRc12=RdRc22,
(34)

where Rc1 is Rc in the study of Nagatsu and De Wit,28 whereas Rc2 is Rc in this study. From Eq. (34), the following relations are obtained:

RdRa=Rb,Rc2=Rc1+RaRb.
(35)

In the study of Nagatsu and De Wit,28,Rb = 2. Then, Ra is set to 1 to satisfy Eq. (35); Rd = 3 and Rc2 = Rc1 − 1 are obtained. Actually, under the condition Rc2 = Rc1 − 1, the same shape of the viscosity profile is obtained as in the study of Nagatsu and De Wit.28 If we set Rb = 3 and Ra = 0.5, Rd = 3.5 and Rc2 = Rc1 − 2.5 are obtained for Eq. (35). To verify this calculation, the underlying viscosity profiles under the condition of Rb = 3, where the viscosity depends on b and c, and of Ra = 0.5, where the viscosity depends on a, c, and d, are shown in Fig. 6. We confirm that the conditions of Rd = 3.5 and Rc2 = Rc1 − 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 Wit28 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.

FIG. 6.

The notation “NR” is shorthand for “nonreactive.” Underlying viscosity profile for (a) various Rc wherein Rb = 3 and the viscosity depends on b and c and (b) Ra = 0.5 wherein Rd = 3.5 and the viscosity depends on a, c, and d.

FIG. 6.

The notation “NR” is shorthand for “nonreactive.” Underlying viscosity profile for (a) various Rc wherein Rb = 3 and the viscosity depends on b and c and (b) Ra = 0.5 wherein Rd = 3.5 and the viscosity depends on a, c, and d.

Close modal

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

2ψ=Rbψxbx+ψyby+by+Rcψxcx+ψycy+cy.
(36)

By substituting Eq. (27) into Eq. (36), we obtain the following equation for the non-reactive case:

2ψ=Rbψxbx+ψyby+by=Rbde2ψx+ψy+1=Rb2de  for all η,
(37)

where we denote the spatial derivative of the error function as de and we use the following equation:

de=deψx+ψy+1.
(38)

By substituting Eq. (28) into Eq. (36), we obtain the following equation for the reactive case:

2ψ=Rc2de,η<0RbRc2de,η0.
(39)

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:

2ψ=RdRa2de, for all η,
(40)
2ψ=Ra+RC2+Rd2de,η<0Rc2+Rd2de,η0.
(41)

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 Rb = Rc provides the equivalent vorticity equation for all η in the non-reactive and reactive cases [Eqs. (37) and (39)].

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 Wit28 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 Wit28 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 CO2 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.

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

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.

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/L3 (M is the dimension of mass), and diffusion coefficients of chemical species A, B, and C, DA,B,C, are constant is written as46 

ρAt+uρA=DA2ρArA,
(A1)
ρBt+uρB=DB2ρBrB,
(A2)
ρCt+uρC=DC2ρC+rC,
(A3)

where ρA, ρB, ρc are the mass concentration of chemical species A, B, and C, whose dimensions are M/L3. In addition, rA and rB are the mass consumption rates of chemical species A and B, which are the reactants, while rC is the mass production rate of chemical species C, the product. The dimensions of rA, rB, and rC are M/(L3t), where t is the dimension of time. When Eqs. (A1)–(A3) are divided by the molecular weight of chemical species A, B, and C, MA, MB, and MC, whose dimensions are M/mole, the following equations are obtained:

at+ua=DA2arAMA,
(A4)
bt+ub=DB2brBMB,
(A5)
ct+uc=DC2c+rCMC.
(A6)

Here, a = ρA/MA, b = ρB/MB, and c = ρC/MC. The reaction terms are now expressed with the dimension of moles/(L3t) as follows:

rAMA=rBMB=rCMC=kab.
(A7)

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

rA+rB=MA+MBkab=MCkab=rC,
(A8)

where MA + MB = MC 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.

1.
S.
Hill
, “
Channelling in packed columns
,”
Chem. Eng. Sci.
1
,
247
(
1952
).
2.
P. G.
Saffman
and
G. I.
Taylor
, “
The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid
,”
Proc. R. Soc. London, Ser. A
245
,
312
(
1958
).
3.
R. L.
Chuoke
,
P.
van Meurs
, and
C.
van der Poel
, “
The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media
,”
Trans. AMIE
216
,
188
(
1959
).
4.
A.
Muggeridge
,
A.
Cockin
,
K.
Webb
,
H.
Frampton
,
I.
Collins
,
T.
Moulds
, and
P.
Salino
, “
Recovery rates, enhanced oil recovery and technological limits
,”
Philos. Trans. R. Soc., A
372
,
20120320
(
2014
).
5.
J. M.
Nordbotten
,
M. A.
Celia
, and
S.
Bachu
, “
Injection and storage of CO2 in deep saline aquifers: Analytical solution for CO2 plume evolution during injection
,”
Transp. Porous Media
58
,
339
(
2005
).
6.
B. S.
Broyles
,
R. A.
Shalliker
,
D. E.
Cherrak
, and
G.
Guiochon
, “
Visualization of viscous fingering in chromatographic columns
,”
J. Chromatogr. A
822
,
173
(
1998
).
7.
G. M.
Homsy
, “
Viscous fingering in porous media
,”
Annu. Rev. Fluid Mech.
19
,
271
(
1987
).
8.
K. V.
McCloud
and
J. V.
Maher
, “
Experimental perturbations to Saffman–Taylor flow
,”
Phys. Rep.
260
,
139
(
1995
).
9.
H. A.
Nasr-El-Din
,
K. C.
Khulbe
,
V.
Hornof
, and
G. H.
Neale
, “
Effects of interfacial reaction on the radial displacement of oil by alkaline solutions
,”
Rev. Inst. Fr. Pet.
45
,
231
(
1990
).
10.
A. R.
White
and
T.
Ward
, “
CO2 sequestration in a radial Hele–Shaw cell via an interfacial chemical reaction
,”
Chaos
22
,
037114
(
2012
).
11.
K. R.
Bhaskar
,
P.
Garik
,
B. S.
Turner
,
J. D.
Bradley
,
R.
Bansil
,
H. E.
Stanley
, and
J. T.
LaMont
, “
Viscous fingering of HCl through gastric mucin
,”
Nature
360
,
458
(
1992
).
12.
A.
De Wit
and
G. M.
Homsy
, “
Nonlinear interactions of chemical reactions and viscous fingering in porous media
,”
Phys. Fluids
11
,
949
(
1999
).
13.
A.
De Wit
and
G. M.
Homsy
, “
Viscous fingering in reaction-diffusion systems
,”
J. Chem. Phys.
110
,
8663
(
1999
).
14.
M.
Jahoda
and
V.
Hornof
, “
Concentration profiles of reactant in a viscous finger formed during the interfacially reactive immiscible displacements in porous media
,”
Powder Technol.
110
,
253
(
2000
).
15.
Y.
Nagatsu
and
T.
Ueda
, “
Effects of reactant concentrations on reactive miscible viscous fingering
,”
AIChE J.
47
,
1711
(
2001
).
16.
J.
Fernandez
and
G. M.
Homsy
, “
Viscous fingering with chemical reaction: Effect of in-situ production of surfactants
,”
J. Fluid Mech.
480
,
267
(
2003
).
17.
Y.
Nagatsu
,
K.
Matsuda
,
Y.
Kato
, and
Y.
Tada
, “
Experimental study on miscible viscous fingering involving viscosity changes induced by variations in chemical species concentrations due to chemical reactions
,”
J. Fluid Mech.
571
,
475
(
2007
).
18.
T.
Podgorski
,
M. C.
Sostarecz
,
S.
Zorman
, and
A.
Belmonte
, “
Fingering instabilities of a reactive micellar interface
,”
Phys. Rev. E
76
,
016202
(
2007
).
19.
Y.
Nagatsu
,
S.-K.
Bae
,
Y.
Kato
, and
Y.
Tada
, “
Miscible viscous fingering with a chemical reaction involving precipitation
,”
Phys. Rev. E
77
,
067302
(
2008
).
20.
Y.
Nagatsu
,
A.
Hayashi
,
M.
Ban
,
Y.
Kato
, and
Y.
Tada
, “
Spiral pattern in a radial displacement involving a reaction producing gel
,”
Phys. Rev. E
78
,
026307
(
2008
).
21.
T.
Gérard
and
A.
De Wit
, “
Miscible viscous fingering induced by a simple A + B → C chemical reaction
,”
Phys. Rev. E
79
,
016308
(
2009
).
22.
Y.
Nagatsu
,
Y.
Kondo
,
Y.
Kato
, and
Y.
Tada
, “
Effects of moderate Damköhler number on miscible viscous fingering involving viscosity decrease due to a chemical reaction
,”
J. Fluid Mech.
625
,
97
(
2009
).
23.
S. H.
Hejazi
and
J.
Azaiez
, “
Non-linear interactions of dynamics reactive interface in porous media
,”
Chem. Eng. Sci.
65
,
938
(
2010
).
24.
Y.
Nagatsu
,
C.
Iguchi
,
K.
Matsuda
,
Y.
Kato
, and
Y.
Tada
, “
Miscible viscous fingering involving viscosity changes of the displacing fluid by chemical reactions
,”
Phys. Fluids
22
,
024101
(
2010
).
25.
S. H.
Hejazi
and
J.
Azaiez
, “
Hydrodynamic instability in the transport of miscible reactive slices through porous media
,”
Phys. Rev. E
81
,
056321
(
2010
).
26.
S. H.
Hejazi
,
P. M. J.
Trevelyan
,
J.
Azaiez
, and
A.
De Wit
, “
Viscous fingering of a miscible reactive A + B → C interface: A linear stability analysis
,”
J. Fluid Mech.
652
,
501
(
2010
).
27.
Y.
Nagatsu
,
Y.
Kondo
,
Y.
Kato
, and
Y.
Tada
, “
Miscible viscous fingering involving viscosity increase by a chemical reaction with moderate Damköhler number
,”
Phys. Fluids
23
,
014109
(
2011
).
28.
Y.
Nagatsu
and
A.
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
).
29.
L. A.
Riolfo
,
Y.
Nagatsu
,
S.
Iwata
,
R.
Maes
,
P. M. J.
Trevelyan
, and
A.
De Wit
, “
Experimental evidence of reaction-driven miscible viscous fingering
,”
Phys. Rev. E
85
,
015304(R)
(
2012
).
30.
H.
Alhumade
and
J.
Azaiez
, “
Stability analysis of reversible reactive flow displacements in porous media
,”
Chem. Eng. Sci.
101
,
46
(
2013
).
31.
Q.
Yuan
and
J.
Azaiez
, “
Cyclic time-dependent reactive flow displacements in porous media
,”
Chem. Eng. Sci.
109
,
136
(
2014
).
32.
P. H.
Bunton
,
M. P.
Tullier
,
E.
Meiburg
, and
J. A.
Pojman
, “
The effect of a crosslinking chemical reaction on pattern formation in viscous fingering of miscible fluids in a Hele–Shaw cell
,”
Chaos
27
,
104614
(
2017
).
33.
Z.
Niroobakhsh
,
M.
Litman
, and
A.
Belmonte
, “
Flow instabilities due to the interfacial formation of surfactant–fatty acid material in a Hele–Shaw cell
,”
Phys. Rev. E
96
,
053102
(
2017
).
34.
S.
Stewart
,
D.
Marin
,
M.
Tullier
,
J.
Pojman
,
E.
Meiburg
, and
P.
Bunton
, “
Stabilization of miscible viscous fingering by a step growth polymerization reaction
,”
Exp. Fluids
59
,
114
(
2018
).
35.
V.
Sharma
,
S.
Pramanik
,
C.-Y.
Chen
, and
M.
Mishra
, “
A numerical study on reaction-induced radial fingering instability
,”
J. Fluid Mech.
862
,
624
(
2019
).
36.
R.
Tsuzuki
,
T.
Ban
,
M.
Fujimura
, and
Y.
Nagatsu
, “
Dual role of surfactant-producing reaction in immiscible viscous fingering evolution
,”
Phys. Fluids
31
,
022102
(
2019
).
37.
D. M.
Escala
,
A.
De Wit
,
J.
Carballido-Landeira
, and
A. P.
Muñuzuri
, “
Viscous fingering induced by a pH-sensitive clock reaction
,”
Langmuir
35
,
4182
(
2019
).
38.
C.
Rana
and
A.
De Wit
, “
Reaction-driven oscillating viscous fingering
,”
Chaos
29
,
043115
(
2019
).
39.
S.
Sin
,
T.
Suekane
,
Y.
Nagatsu
, and
A.
Patmonoaji
, “
Three-dimensional visualization of viscous fingering for non-Newtonian fluids with chemical reactions that change viscosity
,”
Phys. Rev. Fluids
4
,
054502
(
2019
).
40.
R.
Tsuzuki
,
Q.
Li
,
Y.
Nagatsu
, and
C. Y.
Chen
, “
Numerical study of immiscible viscous fingering in chemically reactive Hele–Shaw flows: Production of surfactants
,”
Phys. Rev. Fluids
4
,
104003
(
2019
).
41.
P.
Shukla
and
A.
De Wit
, “
Influence of the Péclet number on reactive viscous fingering
,”
Phys. Rev. Fluids
5
,
014004
(
2020
).
42.
Y.
Nagatsu
, “
Viscous fingering phenomena with chemical reactions
,”
Curr. Phys. Chem.
5
,
52
(
2015
).
43.
A.
De Wit
, “
Chemo-hydrodynamic patterns in porous media
,”
Philos. Trans. R. Soc., A
374
,
20150419
(
2016
).
44.
A.
De Wit
, “
Chemo-hydrodynamic patterns and instabilities
,”
Annu. Rev. Fluid Mech.
52
,
531
(
2020
).
45.
C. T.
Tan
and
G. M.
Homsy
, “
Simulation of nonlinear viscous fingering in miscible displacement
,”
Phys. Fluids
31
,
1330
(
1988
).
46.
R. B.
Bird
,
W. E.
Stewart
, and
E. N.
Lightfoot
,
Transport Phenomena
, 2nd ed. (
John Wiley & Sons
,
2002
).