Pore-scale simulation of miscible viscous fingering with dissolution reaction in porous media

Global climate change is happening but may be mitigated by the technology of geological carbon dioxide (CO 2 ) sequestration. To gain comprehensive insights into this approach, we perform pore-scale simulations of displacement between two miscible ﬂuids in porous media using a new multiple-relaxation-time lattice Boltzmann model. This study marks the ﬁrst attempt to investigate viscous ﬁngering dynamics in miscible displacement, considering the coexistence of viscosity contrast and dissolution reaction. Simulation results capture different ﬁn-gering patterns that depend on dissolution (Damk € ohler number Da), diffusion (Peclet number Pe), and viscosity contrast (viscosity ratio R). From simulations of unstable viscous ﬂows, dissolution is found to delay ﬁngering onset, slow down ﬁngering propagation, and inhibit or reinforce the late-stage ﬁngering intensity. In simulations with stable viscosity contrasts, the displacement features ﬁngering phenomena when dissolution is fast enough. In addition, we conduct a parametric study to assess the impact of Pe, R, and Da. The results suggest that increasing Pe or R destabilizes ﬁngering, but increasing Da ﬁrst suppresses and gradually intensiﬁes ﬁngering. Finally, for every ﬁxed Da, we determine the phase boundary between stable and unstable regimes in a Pe–R phase plane. A uniﬁed scaling law is developed to approximate boundary lines obtained under different Da values. By comparing reactive and nonreactive cases, we classify four distinct regimes: stable, unstable, reactive stable, and reactive unstable. These pore-scale insights are helpful in


I. INTRODUCTION
Global climate change is associated with the growing carbon dioxide (CO 2 ) emission to the atmosphere, thereby directing efforts toward limiting the anthropogenic CO 2 emission. A promising solution is to capture CO 2 from power plants and inject it into geological reservoirs for long-term sequestration. 1 After injection, supercritical CO 2 is less viscous than the displaced brine, and fluid mobility decreases in the flow direction, which induces viscous fingering instability at the interface between the two miscible fluids. 2 Meanwhile, the injected CO 2 acidifies the subsurface brine and dissolves carbonate minerals. Such mineral dissolutions consume CO 2 and modify local porosity (mobility), hence variations in fingering dynamics. 1,3 In this sequestration process, viscous fingering not only promotes the mixing and trapping of CO 2 but also creates potential CO 2 leakage pathways. 4 Therefore, a successful CO 2 sequestration process requires pre-estimations of the displacement stability between the injected CO 2 and the host brine in porous reservoirs. A better understanding of miscible viscous fingering with dissolution reaction in porous media is, thus, of great importance.
Viscous fingering arises when a high-mobility fluid displaces a low-mobility one, forming the finger-shaped penetration of the displacing fluid into the displaced one. 5 This unstable displacement pattern has been extensively analyzed in the existing literature. One of the pioneering works was conducted by Saffman and Taylor who identified fingering phenomena in immiscible displacement. 6 Researchers subsequently extended to study miscible viscous fingering and summarized different stages (i.e., fingering initiation, growth, interaction, and decay) in the full life cycle of viscous fingering. 4 To be connected with industrial applications, miscible fingering dynamics was investigated in the presence of the magnetic field 7 and nanoparticles. 8 The effects of media heterogeneity were studied and showed intensified fingering in high-permeability porous layers (channeling paths). 9,10 In addition to these nonreactive studies, efforts were devoted to modeling fingering dynamics in miscible displacement between two reactive fluids. Numerical results demonstrated chemically driven viscous fingering and new fingering phenomena (droplet and tip splitting). 11,12 Fingering behaviors were also investigated with adsorption reactions on porous matrices, suggesting that adsorption could be a useful strategy to manipulate viscous fingering. 13,14 Recent attention turned to classifying different fingering patterns in parameter spaces, which were determined by relative strengths of convection, diffusion, and chemical reaction. 15,16 These classifications could play guiding roles in practical applications as they revealed possible measures, such as choosing suitable reactants, for controlling viscous fingering. 15 In the above research studies, the fluid mobility contrast, acting as the determination of viscous fingering, comes from differences between fluid viscosities. On the other hand, the mobility gradient can develop in porous media as dissolution reaction occurs on porous matrices. Explicitly, the reaction consumes chemical solutes and dissolves solid matrices, which ultimately increases media porosity (fluid mobility) behind the moving interface. 3 Such mineral dissolutions make a fluid in a high-mobility area displace another fluid in a lowmobility zone, thereby triggering fingering phenomena (also known as infiltration instability). 17,18 This has long been observed in acid dissolution experiments. Fredd and Fogler were one of the first to experimentally identify different fingering channels under various injection flows and dissolution rates. 19 They suggested the existence of optimal conditions that produced channels with the minimum pore volumes for fluid to break through. Dissolution-induced fingering was, then, observed to be enhanced by media heterogeneity and brine pH. 20 In parallel, numerical models were developed to simulate dissolution reaction in porous media. With these simulations, key influencing factors on dissolution-induced fingering were evaluated; 21 parameter diagrams were created to represent fingering patterns; 22 and experimental observations were validated and explained. 23 These existing studies have enriched our knowledge about miscible fingering instability induced by viscosity contrast or dissolution reaction in porous media. However, they were carried out at the macroscopic scale, thus obtaining effective transport and reactive parameters via empirical correlations or volume-averaged techniques. In this way, large uncertainties would be introduced if without pore-scale information. 24 Moreover, the determination of media permeability is complex and not a simple function of media porosity. 25 For these reasons, pore-scale simulations of miscible displacement were performed. A random network model was proposed to describe displacement with chemical dissolution in porous media, with the formation of fingering being captured. 26 The finite difference and the random walk methods were combined to model reactive transport with solid structure evolutions at the pore scale. The results identified different dissolution instabilities based on relative magnitudes of convection, diffusion, and dissolution. 27,28 Nevertheless, these conventional methods discretize macroscopic continuum equations and cannot easily handle discontinuous variables at complex fluid-solid interfaces, limiting their applications for flows in porous media with complicated structures. 29 Over the past three decades, the lattice Boltzmann (LB) method, based on mesoscopic kinetic equations, has been developed and become a powerful alternative to conventional numerical solvers for transport phenomena at the pore scale. This is attributed to its attractive advantages, such as simple implementation, high parallelism, and ability to handle complex physics and boundary conditions. 30,31 Consequently, LB models have been developed for simulating either viscous fingering or dissolution reaction in porous media at the pore scale. On one hand, Liu and Guo 32 developed an LB model to quantify fluid mixing in both homogeneous and heterogeneous porous media, showing viscous fingering introduced disorders in fluid velocity and enhanced fluid mixing. Such a promoting effect of fingering was proved to be suppressed by chemical reactions that could sharp fingering interfaces. 33 For miscible viscous fingering, the temperature and nanoparticles were found to play imperative roles as they modified fluid viscosity. 34 On the other hand, Kang et al. 29 proposed a porescale LB model to study dissolution reactions, with evolutions of the porous structure being tracked by the volume of the pixel scheme. Their results qualitatively verified experimental results in Ref. 19 and suggested the reliability of the LB method for modeling chemical dissolution in porous media at the pore scale. After this first step, they further predicted different dissolution patterns and porosity-permeability relationships under varying flow and dissolution conditions, 35 and extended to model multicomponent systems that were associated with geological CO 2 sequestration. 36 In summary, the above studies have provided detailed insights into fingering dynamics introduced by either viscosity contrast or dissolution reaction in porous media. To the best of our knowledge however, no pore-scale study has been devoted to investigating fingering dynamics considering the coexistence of viscosity contrast and chemical dissolution. In applications such as long-term CO 2 sequestration, these two factors usually take place simultaneously and have significant effects on viscous fingering instability. Although existing simulations have modeled fingering with both viscosity contrast and adsorption reaction, they applied macroscopic assumptions and ignored the dissolution of the porous matrix. 13,14 Therefore, to fill this gap, a multiple-relaxation-time (MRT) LB model is built to simulate miscible displacement in porous media at the pore scale. We focus on studying fingering dynamics in the presence of both viscosity contrast and dissolution reaction.

II. MATHEMATICAL AND PHYSICAL MODELS
In this study, miscible displacement between two fluids (labeled as 1 and 2) is investigated in porous media at the pore scale. As depicted in Fig. 1, a 2D homogeneous structure with length l x and width l y is built. 37 This porous network contains a staggered array of circular and nonreactive grains, which are homogeneously coated by

ARTICLE
scitation.org/journal/phf reactive solid B. Pore spaces between these grains are initially filled with fluid 2 of viscosity l 2 . Another fluid 1 of viscosity l 1 is injected from the left inlet by a driving force F ¼ ðF x ; 0Þ to displace fluid 2. These two fluids are considered miscible, incompressible, isothermal, and neutrally buoyant. During such a displacement, viscosity fingering is expected to occur if the displacing fluid 1 is less viscous than fluid 2 (l 1 < l 2 ). In the meantime, the injected fluid 1 contains a solute A in concentration c ¼ c 1 . Solute A can act as a dissolving component to react with solid B following the A þ B ! P scheme, with the reaction product P being in the dissolved form. Note that, fluid 2 does not contain solute A and its concentration of A is c ¼ c 2 ¼ 0. Thus, in the course of displacement, such a dissolution reaction takes place as solute A from fluid 1 comes into contact with solid B. The reaction rate F r is estimated as 38 where k is the kinetic reaction constant and c I is the concentration of A at interface I between fluid 1 and solid B. c s ¼ c pr =k eq is the equilibrium concentration, with c pr and k eq being the reaction product concentration and the effective equilibrium constant, respectively. In this study, we set c s ¼ 0 for the irreversible dissolution reaction. As solid B is dissolved out gradually, the update of porous geometry is tracked by 38 where V B and V Bm are the volume and molar volume of solid B, respectively, and S B is the reactive surface area. Such a dissolution reaction alters media porosity and permeability, leading to variations in fluid mobility and interface instability. Based on the above introductions, governing equations for fluid flow and concentration transport in pore spaces can be built as where u ¼ ðu; vÞ, p, q 0 , and c are the fluid velocity, pressure, density, and concentration (of species A), respectively. t is the time, and D is the diffusion coefficient. l ¼ q 0 is the fluid dynamic viscosity with being the kinematic viscosity. The effect of mass transfer on fluid viscosity is described by the widely used exponential law as 33 where R ¼ ln ðl 2 =l 1 Þ is the log-viscosity ratio and l 1 and l 2 are the dynamic viscosities of fluid 1 (c ¼ c 1 ) and fluid 2 (c ¼ c 2 ¼ 0), respectively. To complete the present governing equations, dissolution reaction at interface I is described by boundary conditions as 38 No À slip velocity : Species conservation : n Á Drc I ¼ F r ; where n is the interface normal pointing to the fluid phase. For nonreactive cases, Eq. (8) becomes rc I ¼ 0. By introducing the characteristic length L ¼ l y , velocity U ¼ 2 =L, and concentration c 1 , dimensionless parameters marked by asterisks can be derived as Key characteristic numbers of this problem are the viscosity ratio (R), the Reynolds number (Re), the Peclet number (Pe), and the Damk€ ohler number (Da). To solve the above governing equations (3)-(5) and the interface dissolution reaction (7) and (8), the MRT model is used as the basis. 31 This model has been verified to be capable of avoiding the unphysical dependence of permeability on viscosity for pore-scale simulations. Recently, we have proposed a two-dimensional ninevelocity (D2Q9) MRT LB model to simulate coke combustion in porous media. 39 In this study, we extend this model for pore-scale simulations of miscible viscous fingering with dissolution reaction in porous media. Details of the new MRT LB model can be found in Appendix A.

III. RESULTS AND DISCUSSION
We proceed to study viscous fingering dynamics in the context of miscible displacement with dissolution reaction. To conduct pore-scale simulations in the homogeneous medium as built in Fig. 1, we set geometric parameters as l x ¼ 1, l y ¼ 0:497, d ¼ 0.0063, d ¼ 0:001, and r x ¼ r y ¼ 0:015, respectively, leading to the medium porosity as / ¼ 0:508 and the initial volume fraction of solid B as / B;0 ¼ 0:232. In regard to boundary conditions, the top and bottom boundaries are periodic, zero-gradient conditions are applied for all the scalars at the outlet, and the fixed viscosity and concentration of fluid 1 are employed at the inlet. Moreover, the Reynolds number of fluid 2 and the driving force are fixed as Re ¼ 1:0 and F x ¼ 0:09, respectively. Different values of viscosity ratio R, Peclet number Pe, and Damk€ ohler number Da are chosen to change modeling conditions. Before proceeding further, grid convergence tests have been carried out. A mesh of size N x Â N y ¼ 1920 Â 954 is selected to describe the medium in Fig. 1, with other parameters being d ¼ 12d respectively. More details about model validation and grid convergence tests can be found in the supplementary material. Each of the following simulations will be continued until the front position of fluid 1, l f , reaches the position 0:75l x . Note that, throughout the text, we set all the parameters in lattice units and refer to the front position of fluid 1 as the fluid front.

A. General viscous fingering patterns
A range of pore-scale simulations are first performed to investigate general patterns and dynamics of miscible viscous fingering in the presence of both viscosity contrast and dissolution reaction. Our numerical results reveal different fingering types that are determined by the coupled effects of three factors, namely, viscosity contrast (R), fluid diffusion (Pe), and chemical dissolution (Da). As an example, two sets of simulations are discussed to illustrate these fingering patterns, including cases Ia-Ic with R ¼ 3; Pe ¼ 30, and Da ¼ 0; 0:01; 0:3 and cases IIa-IIc with R ¼ À0:5; Pe ¼ 30, and Da ¼ 0; 0:01; 0:3. These three Da values are selected to represent no dissolution, slow dissolution, and fast dissolution, respectively.
We first consider cases Ia-Ic where the viscosity contrast between fluids 1 and 2 is classically understood as unstable (R ¼ 3). For each test, Fig. 2  In spite of this similar development pattern, the dissolution of solid B is different among the three cases. As an illustration, Fig. 3 provides zoom-in views of concentration and residual solid B distributions. In case Ia with no dissolution [ Fig. 3(a)], solid B keeps its initial distribution and media porosity remains unchanged during the displacement. With the introduction of dissolution in cases Ib and Ic   contribute to the delayed viscous fingering appearance. Additionally, as Da increases from cases Ib to Ic, the fingering delay period lasts longer, and thus, the suppressing effect of dissolution is stronger.

Physics of Fluids
On the other hand, at the late stage, the influence of dissolution becomes nonmonotonic. That is, compared with case Ia [Figs. 2(a2)-2(a4)], chemical dissolution tends to either stabilize [Figs. 2(b2)-2(b4)] or destabilize [Figs. 2(c2)-2(c4)] viscous fingering, depending on the dissolution rate (Da). In case Ia, fingers grow and merge with neighboring ones, forming two obvious fingers during the late stage. This is, however, absent in case Ib, where fingers fade away with time. It is attributed to the fact that the slow dissolution keeps smoothening the fluid viscosity contrast and suppressing fingering development. In case Ic, although the onset of fingering is delayed, the fast dissolution of solid B, in turn, brings in a high-mobility region saturated with the displacing fluid 1 behind the moving interface. This subsequently enlarges the fluid mobility contrast and results in the enhanced fingering dynamics (i.e., long fingering and tip splitting). Such a destabilizing effect of dissolution is not observed in case Ic because, as mentioned above, the expansion of the high-mobility area is obviously slower than the propagation of fluid 1. In addition, profiles of transversely averaged concentration c Ã exhibit plateaus in cases Ia and Ic (Fig. 2). These plateaus emerge because fingers contain almost the same quantity of A. By comparing the plateau length, it is also concluded that viscous fingering is reinforced by fast dissolution but inhibited by slow dissolution at the late development stage.
To quantitatively grasp fingering dynamics in cases I, Fig. 5 plots temporal evolutions of front position l f =l x , mixing length l m =l x , volume fraction of fluid 1 Rv 1 ¼ V 1 =V p , and residual solid B R B ¼ / B =/ B;0 . Parameters V 1 , V p , and / B are the volume occupied by fluid 1, the volume of pore spaces, and the volume fraction of B at a given time, respectively. Note that, based on profiles of transversely averaged concentration c Ã , l f is defined as the position with c Ã ¼ 0:01 and l m is calculated as the distance between c Ã ¼ 0:01 and c Ã ¼ 0:99. The curves of l f =l x in Fig. 5(a) show that nonreactive case Ia propagates faster than two reactive cases Ib and Ic. It is because chemical dissolution consumes species A, which contributes to decreasing viscosity (or increasing mobility) of the displacing fluid 1. To understand the behaviors of l m =l x as a function of time in Fig. 5(b), we recall that l m =l x grows in a diffusive tendency during the initial period but starts to depart from this tendency once viscous fingering appears, exhibiting two different growth velocities. 4,33 So, it is observed from profiles of l m =l x that dissolution delays the onset of fingering in reactive cases Ib and Ic and even suppresses the late-stage fingering growth in case Ib. The results of Rv 1 and R B in Figs. 5(c) and 5(d) suggest that the dissolution of B (or the expansion of the high-mobility area) is slower than the propagation of fluid 1 in case Ib. Finally, by comparing final states (with l f =l x ¼ 0:75) of cases Ia-Ic, we can see that the fast dissolution case Ic features the strongest viscous fingering with the longest mixing length (or the fingering length). These results quantitatively confirm the above observations from Fig. 2.
For a comprehensive understanding, we, then, focus on cases IIa-IIc with a stable viscosity gradient R ¼ À0:5. Figure 6 shows the spatial-temporal evolutions of concentration c Ã and transversely averaged concentration c Ã and porosity /. In cases IIa and IIb, diffusion dominates the displacement and fluid interface remains almost undeformed. Once fast dissolution comes into play in case IIc, the fluid interface is distorted with time and finally becomes fingered, regardless of the stable viscosity contrast. As explained in cases I, it is because fast dissolution results in an unstable situation of high-mobility fluid 1 displacing low-mobility fluid 2, hence the triggered viscous fingering (or infiltration instability). Different displacement dynamics between cases IIb and IIc suggest that the onset of fingering requires sufficiently fast dissolution. In addition, compared with cases I, fingering is less intense here and no obvious plateau is observed in profiles of c Ã , indicating the importance of viscosity contrast to viscous fingering. As conducted in cases I, these observations are verified by quantitative profiles as well, which are not provided here for brevity.
In conclusion, for viscous unstable cases (R > 0), we find that dissolution delays the onset of viscous fingering. During the late stage of fingering development, however, dissolution may stabilize or destabilize viscous fingering, depending on the dissolution rate (Da). The suppressed viscous fingering under slow dissolution is a novel finding of our work. On the other hand, in viscous stable situations with fluid viscosity decreasing along the propagation direction (R < 0), fast dissolution is capable of triggering the development of viscous fingering, which, however, is less intense than that in viscous unstable cases.

B. Effects of Pe and R
After discussing general fingering dynamics and patterns, the effects of Pe and R are investigated. Simulations are performed at a fixed Da ¼ 0:3 and a wide range of R and Pe values. Figure 7(a) plots the concentration contour for each test at the time instant when the fluid front moves to l f ¼ 0:6l x . A range of fingering dynamics are visible from these concentration distributions, including no interface instability, marginal fingering, and intense fingering. In tests with small R or Pe, the fluid interface remains undeformed as the stable viscosity contrast or fluid diffusion dominates the system. Nevertheless, with increasing R or Pe, an interface deformation tendency and a transition from planar interface to viscous fingering are noticed. These observations are in qualitative agreement with existing results of miscible viscous fingering without dissolution. 16  Moreover, Fig. 7(b) measures the time period t Ã p required for the fluid front propagating to the position l f ¼ 0:6l x . The smaller t Ã p value indicates the faster displacing speed. All cases are observed to share a common tendency that, with increasing R, the displacing speed increases monotonically. It is because the larger R value represents the less viscous fluid 1. In contrast, with ascending Pe, the displacing speed is found to decrease at first and turns to increase at later times. It can be explained in terms of the stability condition of the fluid interface: with the increase in Pe at each fixed R, the fluid interface transits from stable diffusion to viscous fingering. In the diffusion regime, larger Pe denotes slower diffusion and thus the decreased displacing speed. Once viscous effects take over, fingering development and fluid front propagation are accelerated by ascending Pe. It indicates that the maximum t Ã p at each fixed R locates near the transition point from stable diffusion to unstable fingering.
To quantify the effects of R and Pe, we provide temporal evolutions of front position l f =l x , mixing length l m =l x , volume fraction of fluid 1 Rv 1 , and residual solid B R B , for tests with fixed Da   Figs. 9(a) and 9(b). This is because, as shown in Fig. 7(a), viscous fingering is intensified as Pe or R ascends, causing less volume being occupied by fluid 1 and more solid B being left undissolved behind the fluid front. On the other hand, the curves of l m;e =l x generally ascend with increasing Pe or R, suggesting the intensified viscous fingering. A different tendency should be noted: the profile of l m;e =l x decreases with Pe when Pe is

ARTICLE
scitation.org/journal/phf relatively small. It is because diffusion dominates the displacement and larger Pe represents slower fluid diffusion and smaller l m;e =l x . This declining tendency, however, lasts a small range of Pe values. Then, l m;e =l x turns to increase with Pe, indicating that the system transits into the unstable state with viscous fingering. Therefore, a local minimum in l m;e =l x is noticed, and the corresponding Pe can roughly represent the threshold for the unstable regime.
In general, we have both qualitatively and quantitatively demonstrated that fingering intensity boosts with increasing R or Pe. To grasp the distribution of system stability conditions in a phase plane spanned by Pe and R, we further compare the profiles of l m for fixed Pe and Da but varying R values, as conducted in Fig. 8(f). To stand for the stable state, a reference curve of mixing length l m;r is measured vs time at R ¼ À2:5 for each pair of Pe and Da. The viscosity ratio R ¼ À2:5 is selected after checking that it can make the system stay stable despite the values of Pe and Da. Consequently, we compared l m with l m;r to classify the parameter set (Pe; Da; R) as stable or unstable. If l m coincides with l m;r even up to the final time (when l f ¼ 0:75l x ), the corresponding case with (Pe; Da; R) is identified as stable; otherwise, it is unstable. Accordingly, for Da ¼ 0:3, the obtained stable and unstable cases are summarized in the Pe-R space in Fig. 9(c). The results show that, along constant Pe (or R) lines, the system undergoes a transition from stable planar interface to unstable viscous fingering. This reveals the existence of a critical Pe c at a given R (or a critical R c at a fixed Pe) for the occurrence of viscous fingering. Thus, a boundary line, on which critical parameter pairs (Pe c ; R c ) are distributed, is identified to separate the Pe-R space into a stable and an unstable regimes.

C. Effects of Da
We further explore the effects of the dissolution rate (Da) on fingering dynamics. Considering typical fingering phenomena with varying Da values, which are discussed in Sec. III A, we focus here on quantitative results. For cases with Pe ¼ as Da increases, either l m;e =l x or Rv 1;e gradually approaches an ultimate value, and then, the destabilizing effect of dissolution is reflected by tip splitting phenomena. On the other hand, the results of R B;e demonstrate, as expected, that a larger dissolution rate (Da) brings in lesser residual B. In addition, Fig. 10 compares the time period t Ã e required for the fluid front moving to l f ¼ 0:75l x in each case, showing that t Ã e accelerates with Da. As explained above, this is because dissolution consumes species A and slows down the development speed of fingering. Again, these findings quantitatively confirm observations in Sec. III A. Similar simulations and analyses are conducted for varying Da values at other pairs of (Pe; R), and the results in Fig. 10 remain almost unaffected in terms of qualitative tendency.
We, then, conduct the analysis mentioned in Sec. III B to summarize the obtained stable and unstable cases in the Pe-R parameter space for different values of Da. To grasp the distribution of typical fingering patterns, boundary lines for Da ¼ 0:0; 0:01; 0:3 are first provided in Fig. 11. Similar to results in Fig. 9(c), each boundary line divides the space into a stable area and an unstable one. In addition, the curve for nonreactive cases (Da ¼ 0:0) is sandwiched between the two reactive profiles (Da ¼ 0:01; 0:3). Under these three lines, four different stability regimes are further obtained: unstable (I), stable (II), reactive stable (III), and reactive unstable (IV). As an example, truncated concentration contours are provided in Fig. 11 to help visualize fingering patterns in these four regions. It is emphasized that the reactive stable regime, as a novel finding, stands for stable displacements even under nonreactive unstable conditions. These results help to enrich our understanding about the effects of dissolution on viscous fingering. Furthermore, Fig. 12(a) shows boundary lines between stable and unstable regions for a wide range of Da values in the Pe-R parameter space. It is observed that, with increasing Da, the unstable region spans over a smaller range of both Pe and R than the nonreactive case at first but turns to span a larger area after a threshold value (approximately at Da c ¼ 0:01). In addition, regime boundary lines for different Da values follow a unified scaling relation as R c a þ b $ Pe À0:6 c , with correlations a and b being related to the dissolution rate (Da). Using this scaling relation, we may predict the stability of miscible displacement

IV. CONCLUSIONS
Based on a new multiple-relaxation-time reactive lattice Boltzmann model, pore-scale simulations of miscible displacement between two fluids have been performed in a homogeneous porous medium. Considering that finger-shaped propagation affects the storage efficiency of the displacing fluid, this study has investigated viscous fingering dynamics in situations with both viscosity contrast and dissolution reaction. Specifically, a chemical component is introduced in the displacing fluid to dissolve porous matrices and cause the evolution of the porous structure. Simulations have been conducted over a wide range of values in the Damk€ ohler number Da, Peclet number Pe, and viscosity ratio R. The obtained results demonstrate different fingering patterns. In situations of a less viscous fluid displacing a more viscous one (R > 0), dissolution decelerates both the onset and the development of viscous fingering. During the late development stage, viscous fingering gradually fades away in slow-dissolution cases but becomes intensified and features long fingers and tip splitting under fast-dissolution conditions. For cases with stable viscosity contrast (R < 0), fast dissolution is observed to increase media porosity behind the moving interface and thus trigger viscous fingering. In addition, the effects of Pe; R, and Da have been quantitatively analyzed. The results suggest that, as Pe or R increases, more destabilized viscous fingering is evidenced by longer and clear fingers. In contrast, with increasing Da, the fingering intensity is first suppressed and then enhanced by chemical dissolution. By measuring the mixing length, we identified and summarized stable and unstable cases in a parameter space spanned by Pe and R. For each fixed Da, a boundary line is determined in the Pe-R space between a stable and an unstable region. All the boundary lines for different Da values show a similar pattern and can be scaled by a unified empirical law. Furthermore, the boundary line for nonreactive cases (Da ¼ 0) is sandwiched by reactive lines (Da > 0). Thus, four distinct stability regimes are identified as unstable, stable, reactive stable, and reactive unstable.
Based on the present pore-scale simulations, improved physical insights into miscible viscous fingering with dissolution reaction are obtained. The findings have implications for many applications. For example, in the context of geological sequestration of CO 2 or enhanced oil recovery, different storage sites or oil reservoirs should be compared in terms of their chemical components. Sites with slow dissolution reaction should be chosen to help suppress fingering development and thus reduce the CO 2 leakage risk and improve the oil recovery efficiency. Considering that solid matrices and reactive solid B usually distribute randomly in practical applications, further attention will be paid to the influence of porous media heterogeneity on the development of viscous fingering with dissolution reaction.

SUPPLEMENTARY MATERIAL
See the supplementary material for the complete model validation and grid convergence tests.  Here, q p is a variable related to the fluid pressure as p ¼ c 2 2 q p , with c s ¼ e= ffiffi ffi 3 p being the lattice sound velocity. e ¼ d x =d t is the lattice speed and set as e ¼ 1 in this work, with d x and d t denoting the lattice spacing and the time step, respectively. To avoid discrete lattice effects, the distribution function for the driving force F is 31 F For the present D2Q9 model, the transformation matrix M is given as This matrix maps the distribution functions from the physical space w to the moment space asŵ ¼ M Á w. With such a transformation, the evolution equations (A1) and (A2) are implemented in the moment space as