Pore-scale study of miscible density instability with viscosity contrast in porous media

The transport of miscible fluids in porous media is a prevalent phenomenon that occurs in various natural and industrial contexts. However, this fundamental phenomenon is usually coupled with interface instabilities (e.g., viscous/density fingering), which has yet to be thoroughly investigated. In this paper, a multiple-relaxation-time lattice Boltzmann method is applied to study the displacement between two miscible fluids in porous media at the pore scale, with the coexistence of density difference (Rayleigh number Ra), viscosity contrast (R), and injection velocity (U top ). A parametric study is conducted to evaluate the impact of Ra, R, and U top on the flow stability. For a fixed Ra that can trigger density fingering, the increase in R or U top is found to suppress density fingering. Consequently, under a large U top and a moderate R, the density fingering is fully stabilized and the flow follows a stabile pattern. Furthermore, as both R and U top grow to a sufficiently high level, they can jointly trigger viscous fingering. In addition, the increasing Ra shows an enhancing effect on both density fingering and viscous fingering. Finally, by quantitatively analyzing the fingering length ( l m ) and the fingering propagation time ( t e ), five different flow patterns are classified as viscosity-suppressed (I), viscosity-enhanced (II), viscosity-unstable (III), displacement-suppressed (IV), and stable (V) regimes. In a three-dimensional parameter space spanned by Ra, R, and U top , the parameter ranges of the five regimes are determined according to l m and t e . These findings hold a significant value in providing guidance for controlling the flow stability by selecting appropriate operating conditions


I. INTRODUCTION
Carbon dioxide (CO 2 ) emission has been an increasing concern due to its role as a main source of greenhouse gases, which cause global warming, iceberg melting, and sea level increasing. 1,2Among the known technologies for reducing CO 2 emissions, carbon capture and storage (CCS) is one of the most promising methods in terms of cost and scale. 3Deep saline aquifers are the most popular storage sites due to their large space and availability. 4In CCS, CO 2 is injected into the aquifer and gradually dissolved into the brine. 5The CO 2 dissolution slightly increases the brine density, thus introducing a buoyantly unstable stratification of a denser CO 2 -enriched brine on top of the fresh brine. 6,7Such a stratification finally triggers density instability accompanied by strong convection, which plays a significant role in enhancing the dissolution and spreading of CO 2 into the aquifers. 8,9eanwhile, aside from the increasing effect on density, the dissolution of CO 2 also modifies the fluid viscosity. 10The viscosity contrast can greatly affect the onset and development of density fingering. 11herefore, it is essential to study density instability with viscosity contrast in porous media to enhance the understanding of CCS.
3][14][15][16][17][18][19][20][21] Experiments were performed in Hele-Shaw cells, micromodels, sand packs, or cores to identify density instability patterns and analyze the underlying mechanisms. 16,17,22The key factors for accelerating fluid transport and affecting fluid mixing were discussed. 13owever, these experiments were hard to exactly control the configuration of porous media or concentration perturbations, which might greatly influence the instability pattern. 21Moreover, experiments were generally time-consuming, costly, and highly subjective to human and equipment capacities. 13Theoretical analysis played an important role in predicting possible instabilities.4][25][26][27][28][29][30][31] These theoretical analyses focused on predicting the onset of fingering via analyzing the growth rate of initial perturbations, which were able to cover the possible instability scenarios at an initial stage before nonlinear behaviors happen. 189][30][31] However, theoretical studies could not describe the whole-life behaviors of density instability.
Numerical simulations are helpful to understanding the CO 2 dissolution and transport behaviors during CCS.For such a process, numerical studies were conducted to investigate mechanisms of both density and viscous instabilities.For density fingering, the whole-life behaviors of fingering development in porous media were simulated and it was revealed that the nonlinear interactions of fingers exert a substantial influence on both the speed of front propagation and the overall mixing. 18Four development stages were identified as diffusion, early convection, late convection, and convection shutdown. 32The flux growth regime was further divided into two sub-regimes: the early flux growth regime caused by tiny fingers and the late flux growth regime caused by the nonlinear perturbations. 21The nonlinear interactions between chemical reactions and density fingering were studied. 33Moreover, the effects of perturbations and porous media ignored in experiments were also reported in numerical simulations. 15,24In CCS, the density instability caused by the dissolution of CO 2 results in an acceleration of mass transport between fluids, thus reducing the time required for CO 2 storage. 34As for viscous fingering, the viscosity contrast was found to introduce mobility differences and thus yield viscous fingering phenomena between two displacement fluids. 35The nonlinear phenomena of shielding and spreading were observed in simulations, together with the tip splitting in miscible displacement first reported by Tan and Homsy. 36Furthermore, the calculated metric of mixing length, serving as an indicator of fingering development, was found to grow linearly over time. 379][40][41][42] Larger viscosity ratios were found to delay the fingering onset but result in an increased number of fingers with smaller characteristic wavelength. 43Similar to density fingering, the development of viscous fingering was divided into different stages according to the fingering's full-lifecycle. 42oreover, within the context of CCS, viscous fingering was observed to introduce flow channels for CO 2 leakage, yet simultaneously induce strong convection for enhancing fluid dissolution. 420][41] Chemical reactions were found to have complex effects on viscous fingering by changing both the fluid concentration and the media porosity. 20,38,44,45he above studies provide much insight into density and viscous instabilities in porous media.However, there remain two main problems.On the one hand, all the above research is considered either density fingering or viscous fingering, which included density or viscosity contrasts.In applications of CCS, however, the dissolution of CO 2 modifies both fluid density and viscosity, thus introducing differences in both density and viscosity between the displacing and displaced fluids.There exist some studies on density-driven convection coupled with viscosity contrast through experimental observations. 11,46It was reported that viscosity contrast played an important role in convective mixing by changing the dissolution rate and onset time. 46Different fingering structures and dissolution behaviors were observed for different viscosity ratios.However, the inherent limitations of experiments restricted the available parameters in experiments, and thus the mechanism underlying the density fingering with viscosity contrast was not fully investigated.On the other hand, due to the complicated structure of porous media, current numerical studies were mainly based on volume-averaged scales, like representative elementary volume scale or Darcy scale.A pioneering study was conducted, utilizing a Hartleytransform-based pseudospectral method to simulate the situations when viscosity and density contrasts have opposite effects on stability. 27However, these volume-average studies were likely to introduce errors if without pore-scale information. 47Therefore, a pore-scale study is needed to show the details of density fingering with viscosity contrast.
This study is aimed at providing a pore-scale simulation of density instability with the coexistence of viscosity contrast, thereby getting a deeper understanding of how the viscosity contrast affects instability development and then the fluid dissolution.In addition, the effects of injection velocity of the top layer are also investigated in this paper as an innovative attempt.During the past decades, the lattice Boltzmann method (LBM) has become a powerful solver for simulating fluid flows in porous media at the pore scale.This is attributed to its advantages in parallel computing and handling complex boundaries. 48In this study, a multiple-relaxation-time (MRT) lattice Boltzmann (LB) method is utilized to solve the governing equations and the nonequilibrium extrapolation method, and the halfway bounce-back method are used to deal with the boundary conditions.

II. CONFIGURATION
As shown in Fig. 1, the computational domain is a homogeneous 2 D system with length l x and width l y .A set of uniformly distributed solid grains is utilized inside the domain to simulate the porous structure in reality.The detailed location and size settings for these grains are shown in Fig. 1.In order to study the fingering phenomenon under the effects of gravity, viscosity difference, and injection velocity, two fluids (labeled as fluid 1 and fluid 2) are introduced in this system.These two fluids are considered miscible, nonreactive, isothermal, and incompressible.Initially, fluid 2 without solute A suffuses the whole domain.Then, fluid 1 with solute A at a concentration of C 0 is injected into the system from the top.The injection velocity is U 0 .Notably, a concentration disturbance is applied to help trigger the fingering phenomenon.In this research, the dissolution of species A is assumed to modify both the fluid density and viscosity.

III. GOVERNING EQUATIONS
The governing equations for the fluid flow and species transport are described by the Navier-Stokes equations and the convection-diffusion equation as follows: where u ¼ u; v ð Þ is the velocity, q 0 is the density, t is the time, p is the pressure, l ¼ q 0 is the fluid dynamic viscosity with being the fluid kinematic viscosity, C is the concentration, and D is the diffusion coefficient.F is the force applied to the flow, which is introduced by gravity difference.According to the Boussinesq approximation, fluid density is assumed as a constant q 0 except in the force term F. The gravity acceleration term caused by the concentration C is defined as 49 where g is the gravity and b is the concentration expansion coefficient.Moreover, the fluid viscosity l is also decided by the concentration evolution.The viscosity is assumed to follow an exponential law with the local concentration C as 50 where R ¼ lnðl 1 =l 0 Þ is the viscosity ratio, and l 0 and l 1 are the dynamic viscosities of fluid 1 (C ¼ C 0 ¼ 1) and fluid 2 (C ¼ C 1 ¼ 0), respectively.As for the boundary conditions, no-slip and no-flux conditions are applied to the grains inside the system.The lateral boundaries of the system are set as periodic.The top boundary is defined by the prescribed velocity and concentration.The bottom boundary is defined by the no-slip condition and given concentration.In summary, these boundary conditions are described as where x g ; y g ð Þ is the surface of solid grains and U 0 is the injection velocity of fluid 1.
In order to introduce dimensionless parameters, the characteristic length, velocity, and concentration are selected, respectively, as The nondimensional parameters are as The key parameters in this project are the Rayleigh number Ra, the viscosity ratio R, the top flow velocity U top , and the Schmidt number Sc.The nondimensional injection velocity of the top flow U top is called the injection velocity for convenience in Secs.IV A-IV C.
In this work, the multiple-relaxation-time (MRT) LB method is utilized to solve the above governing equations and boundary conditions. 51The details of the models are shown in the Appendix.

IV. RESULTS AND DISCUSSION
In this part, the density instability with viscosity contrast is simulated for different injection velocities in a homogeneous porous medium.For the computational domain in Fig. 1, the geometric parameters are set as l x ¼ 1; l y ¼ 0:672, d ¼ 0:008 06; r x ¼ r y ¼ 0:0181, and the medium porosity is calculated as u ¼ 1 Àpd 2 =2r x r y ¼ 0:670.It is worth mentioning that the present porescale study is able to capture representative features, and all typical phenomena of fingering are observed, including merging, shielding, tip splitting, and finger reinitiation.Therefore, the simulation results at pore scales are of relevance at larger scales.In this study, the Schmidt number is fixed as S c ¼ 100, and different values of Rayleigh number Ra, viscosity ratio R, and injection velocity U top are selected to change the simulation conditions.After grid convergence tests, a mesh of size N x Â N y ¼ 1488 Â 1000 is utilized in this research.It is noted that all the simulations use the same initial concentration distribution to exclude the effects of perturbations. 24Throughout this paper, all the parameters are set in lattice units.Each test continues until the fingering front (i.e., the front position of fluid 1, l f Þ reaches the position l f ¼ 0:75l y .More details about model validation can be found in the supplementary material.

A. General pattern analysis
To get a comprehensive understanding of the density fingering phenomena with viscosity contrast at the pore scale, simulations of fingering patterns are conducted for different parameters.For illustration, two sets of simulation cases are discussed: cases I(a)-I(c) without injection velocity and cases II(a)-II(c) with injection velocity.The detailed simulation parameters are listed in Table I.

Without injection velocity (U top 50)
In cases I(a)-I(c), the injection velocity is excluded from consideration.The Rayleigh number is set as Ra ¼ 10 9 , which is large enough to trigger the density fingering phenomenon.The viscosity ratio R is set as 0, 3, and À3 for cases I(a), I(b), and I(c), respectively.For each case, Fig. 2 shows the contours of density fields q Ã at different time instants t Ã , with each row of contours reaching the same front position l f .
It is found that all three cases display a classical fingering development pattern, regardless of the viscosity contrast value. 21,32,42nitially, the spreading of solute A is dominated by diffusion as a planar.With the accumulation of A in the top layer, the density difference between the A-accumulated top layer and the host fresh layer eventually gets large enough to trigger the onset of tiny fingers.These fingers  keep growing and start to interact nonlinearly with each other.Finally, after experiencing a merging state, several strong fingers are left to dominate the flow.Typical phenomena of fingering in miscible displacement are observed in Fig. 2, like dominant fingers' spreading and shielding and tip splitting. 36lthough these three cases share a similar fingering development pattern, three differences among them are observed due to varied viscosity ratios.First, finger re-initiation and tip splitting are enhanced in case I(c), while they are suppressed in case I(b), which are marked by rectangles and ellipses in Fig. 2, respectively.This stems from the fact that, as R increases, the displaced fluid (i.e., fluid 2) manifests larger viscosity and thus suppresses finger re-initiation and tip splitting.Second, the time period each case takes to reach the same front position l f varies.Comparing the non-dimensional time t Ã for each front position, it is evident that the ascending R slows down the propagation of solute A along the direction of gravity.For example, from cases I(c) via I(a) to I(b), it takes t Ã ¼ 71; 357 and 834, respectively, for the fingering front to reach the position l f ¼ 0:75l y .This is straightforward as a larger R represents a more viscous fluid 2, which in turn causes a larger resistance force as fluid 1 tries to spread into fluid 2. Finally, with the increase in R, the enhanced accumulation of A in the top layer is detected, which subsequently brings about the amplified density jump Dq Ã and buoyance force F.
After qualitative discussions of similarities and differences between the three cases I(a)-I(c), quantitative analyses are further set out to provide a deep understanding of density instability with viscosity contrast.For such a purpose, the horizontally averaged density is calculated as In Fig. 2, the horizontally averaged density q Ã is plotted on the right side of each contour field.From profiles of q Ã in each case, the line slope [e.g., marked as d y q Ã at t Ã ¼ 238 in case I(a)] is found to decrease with time.This is attributed to the fact that fingering spreads downward and accelerates the dissolution of species A. In the meantime, as q Ã starts to decrease along the y direction (or d y q Ã is no longer infinite), an inflection point is observed in each profile of q Ã [e.g., marked by a triangle at t Ã ¼ 24 in case I(a)].Accordingly, the density jump Dq Ã between the A-accumulated top layer and the host fresh layer can be determined from profiles of q Ã .As denoted in Fig. 2, Dq Ã is calculated as the density difference from the inflection point to the host fluid density level.It is apparent that Dq Ã determines the buoyance force F and subsequently the onset and development of density fingering.From Fig. 2, it is confirmed that the rise in R causes a growth in Dq Ã or F. This can be explained from the aspect of force balance.With the growing R, the enlarged viscosity force requires a large Dq Ã or F to be balanced.The relationship between R and Dq Ã (or F) can also be understood from the perspective of Taylor dispersion.Taylor dispersion refers to the phenomenon of effective mixing and diffusion that occurs in a shear flow with a concentration gradient. 52In case I, increasing R hinders the spreading velocity of solute A, which, in turn, weakens the intensity of Taylor dispersion.The weakened dispersion causes larger concentration inside the fingers.Therefore, increasing R generates a large Dq Ã .However, it is notable that the simulation in this work ends at l f =l y FIG. 2. Contours of density q Ã and transversely averaged density q Ã for cases I: (a) R ¼ 0, (b) R ¼ 3; and (c) R ¼ À3.
¼ 0:75 with flow dominated by advection, while from a previous study, 52 Taylor dispersion becomes independent of R at the very late stage when the flow is dominated by dispersion instead of advection.
On the other hand, two important metrics are introduced based on profiles of q Ã to quantify the development of density fingering.
First, the fingering front l f is defined as the position where q Ã ¼ 0:01.Second, the mixing length (or fingering length) l m is calculated as the distance between q Ã ¼ 0:99 and q Ã ¼ 0:01.Temporal evolutions of the front position l f =l y and mixing length l m =l y for cases I(a)-I(c) are depicted in Fig. 3.Each simulation stops when l f =l y ¼ 0:75, and the corresponding time is termed as the end time, t e .In cases I(a)-I(c) with U top ¼ 0, t e is observed to increase with the growing R, which identifies the slow spreading of species A. This is expected since the large R introduces the strong viscous resistance to suppress the fingering development driven by the buoyance force F. As for temporal evolutions of the mixing length l m =l y , profiles are similar to those of l f =l y .The reason is that, in cases I(a)-I(c) with U top ¼ 0, no additional force except for the buoyance force is applied to drive the movement of fluid 1.Therefore, density fingering starts almost from the top boundary and l m is the same as l f .Moreover, it is observed that, after the initial diffusion stage, l m changes to grow at an almost constant rate due to the appearance of density fingering, which is consistent with previous studies. 36,37Accordingly, as denoted in Fig. 3, the slope of l m during the fingering development stage is defined as the fingering growth velocity U m .From profiles of l m , U m is found to have a negative correlation relationship with R in cases I(a)-I(c).This is attributed to the suppressing effect of large R on the fingering development.In cases II(a) and II(b), an additional injection velocity is introduced as U top ¼ 10 À2 , while the Rayleigh number and the viscosity ratio are identical to those of cases I(a)-I(c), namely, Ra ¼ 10 9 and R ¼ 0; 3; À3.Similarly, the contours of density fields q Ã for cases II are displayed in Fig. 4. The fingering development patterns in cases II exhibit a distinct divergence from those observed in cases I.In case II(a) where R ¼ 0, the system maintains stable throughout the simulation life span.Conversely, in cases II(b) and II(c) with non-zero R, the onset and development of density fingering are observed.In case II(b) with a high viscosity ratio R ¼ 3, fingers appear at an early stage.The difference between cases II(a) and II(b) reveals that fingers in case II(b) are introduced by the combination of injection velocity and viscosity ratio, which are thus referred to as viscous fingering.On the other hand, in case II(c) with a small viscosity ratio R ¼ À3, the system remains stable for a long period.Upon the fingering front reaching 0.75 l y , the flow in case II(c) undergoes instability and produces minuscule fingers, as denoted by dotted circles in Fig. 4. Therefore, upon examining the contours of R ¼ À3, 0, and 3, the flow field transitions from an unstable regime to a stable one and ultimately reverts to an unstable state.By comparing with cases I, the flow states in cases II indicate that the injection velocity serves to suppress the development of density fingering but enhance that of viscous fingering as R > 0.
In addition to fingering patterns, the injection velocity is observed to accelerate the spreading of solute A. Upon examining the time instant t Ã when l f reaches the same position, it is apparent that t Ã is significantly smaller in cases II as compared to that in cases I.For example, it takes t Ã ¼ 71 in case I(a) but t Ã ¼ 50 in case II(a) for the fingering front to reach l f ¼ 0:75l y .In addition, a comparison of the fluid transport speed among cases II(a)-II(c) demonstrates that the increasing R decelerates the fluid propagation velocity as R < 0, which is consistent with cases I.In contrast to cases I, as R increases to be R > 0, the higher R changes to accelerate fluid transport and the case II(b) with R ¼ 3 shows the fastest transport velocity among cases II.This discrepancy is attributed to the fact that the injection velocity and large viscosity R ¼ 3 in case II(b) bring about viscous fingering.Fingering is known to expedite fluid transport, and therefore, even though a higher viscosity ratio generates the stronger viscous resistance in case II(b), the intensified spreading by fingering surpasses the delaying effect induced by the viscous resistance.
Following qualitative discussions of fingering properties in cases II, quantitative analyses are conducted to comprehend the understanding of density instability with injection velocity.Again, the horizontally averaged density q Ã is graphed on the right side of each contour plot.For each profile of q Ã in Fig. 4, an inflection point is identified and then the density jumps Dq Ã (or the buoyance force F) can be obtained.Different from cases I, values of Dq Ã in cases II(a)-II(c) are almost the same and reach the maximum density jump Dq Ã m ¼ 1.Although under Dq Ã m , density fingering is weak or even disappeared in cases II(c) and II(a).This verifies that the injection velocity U top impedes the effects of F and subsequently suppresses the density fingering development.
Temporal evolutions of the aforementioned metrics l f and l m are plotted in Fig. 5 to quantify the fingering development with injection velocity in cases II.For each test, the simulation ends at t e with l f ¼ 0:75l y .Profiles of l f are observed to increase with time following a similar trend as those in cases I.Each l f increases progressively until the end time t e .In cases II with U top ¼ 10 À2 , as R increases, t e is found to increase when R < 0 but switches to decreases when R > 0. Different from cases I with the positive correlative between R and t e , this nonmonotonical relationship between R and t e in cases II evidences the enhancing effects of viscous fingering on fluid propagation velocity in case II (b) with R ¼ 3. Another distinction from cases I is that, driven by the additional U top , results for l m exhibit a significant difference from those for l f .In case II(a) (R ¼ 0), l m follows a growth pattern dominated by diffusion with no fingering observed.Differently, l m for case II(b) (R ¼ 3) departs from the diffusive tendency at an early stage and turns to increase rapidly due to the onset of viscous fingering, whereas l m for case II(c) (R ¼ À3) does not deviate from the diffusive trend until a very late stage, indicating the weak density fingering, which is also in good agreement with previous study. 27Therefore, the largest growth velocity of fingering, U m , among cases II(a)-II(c) is found in case II(b) with R ¼ 3.This is attributed to the fact that the injection velocity impedes the onset of density fingering but helps to trigger viscous fingering when R > 0.

B. Effects of R and U top
Following a preliminary analysis of density fingering properties with viscosity contrast and injection velocity, effects of the two key parameters (i.e., R and U top ) are further investigated.To this end, a series of simulations are performed with a fixed Ra ¼ 10 9 and varying R and U top .The simulated density contours at the time instant t e when the fluid front reaches l f ¼ 0:75l y are plotted in Fig. 6(a).These contours reveal a diverse array of fingering dynamics, including enhanced fingering, weak fingering, and stable interface.In tests with a small U top (e.g., 10 À4 ), the fluid remains unstable for various values of R and the increasing R suppresses the fingering intensity.However, as U top increases to be efficiently large (e.g., 10 À2 ), the fluid displacement becomes the dominant mechanism, which suppresses the development of density fingering.Under this condition, the growing R within the range R < 0 aids in restraining the formation of density fingering and even yields a stable fluid interface.For instance, in tests with U top ¼ 10 À2 and R ¼ 0, À1, the system maintains a stable planar fluid interface.This phenomenon can be explained by the opposite effects of viscosity and density contrasts on stability by choosing a favorable injection velocity. 27In the large-U top scenario, the density contrast serves to enhance instability, while the large U top amplifies viscous forces and suppresses density instability.Upon U top reaching a sufficiently large magnitude, the flow can be fully stabilized, and the bound velocity is defined as the critical velocity, U c .For example, in test with R ¼ À1, U c is located between 10 À3 and 10 À2 , since the flow is fully stabilized when U top increases from 10 À3 to 10 À2 .It is notable that U c varies with Ra (density contrast), R (viscosity contrast), and time.
However, as R increases to be positive (i.e., R > 0), the combination of large U top and R is observed to introduce viscous fingering, as displayed in test with U top ¼ 10 À2 and R ¼ 3. On the other hand, at a fixed R, the increasing U top is found to hinder density fingering while triggering viscous fingering if R is efficiently large (i.e., R ¼ 3).
Additionally, Fig. 6(b) displays the time t e required for the fluid front to reach the position l f ¼ 0:75l y .A smaller value of t e signifies a faster spreading speed.It is demonstrated that the rise in U top leads to a reduction in t e , thereby hastening the spreading speed of solute A. This can be explained by the fact that U top accelerates the movement of the displacing fluid 1.In contrast, there is a negative correlation between R and the spreading speed of A when U top is relatively small (i.e., U top 10 À4 ).This negative correlation becomes marginal and even disappears as U top increases (i.e., U top ! 10 À3 ).For instance, with the increasing R at U top ¼ 10 À3 , t e increases at first and reaches the largest value at R ¼ 0 and turns to decrease after this maximum point.This is attributed to the appearance of viscous fingering in tests with large U top and R.
To quantify the effects of R, the mixing length l m and the growth velocity of fingering U m at t e when the front positions reach l f ¼ 0:75 l y are plotted in Fig. 7 as a function of R for different U top .In Fig. 7, the data for U top ¼ 10 À4 is omitted to facilitate readability, given its close resemblance to that for U top ¼ 0. The results reveal that effects of R on fingering properties depend on the magnitude of U top .For a small U top (U top 10 À3 ), a negative correlation between l m and R is evident.However, as U top increases to U top ¼ 10 À2 , l m exhibits a weak negative association with R as R < 0, but a strong positive correlation with R as R > 0. This stems from the fact that, although the large U top or R individually can suppress density fingering, the combination of these two factors can serve as the driving force for the initiation of viscous fingering.With regard to U m , as U top is small (U top 10 À3 ), a negative correlation between U m and R is detected from Fig. 7.This correlation arises from the fact that the large R tends to postpone the onset of fingering.However, results at U top ¼ 10 À2 exhibit a nonlinear association between U m and R. Specifically, U m decreases with the growing R as R < 0, but U m experiences an increase as R > 0. For example, the significant rise in U m was detected at R ¼ 3.This nonlinear correlation further corroborates our earlier analysis that the large R suppresses density fingering, but the large U top and R can trigger the appearance of viscous fingering.
Moreover, to quantify the effects of U top , the mixing length l m and the growth velocity of fingering U m upon a front position of l f ¼ 0:75 l y are plotted in Fig. 8 as a function of U top for different R. The results indicate that an increase in U top generally leads to a reduction in l m , which is attributed to the suppressing effect of U top on density fingering.The results reveal a negative correlation between U top and l m , which further corroborates the previous conclusion that U top generally suppresses fingering.On the other hand, the results reveal that U top has minimal impact on U m , except in tests involving large U top and R. For instance, a considerable surge in U m is detected when U top reaches U top ¼ 10 À2 for R ¼ 3.This phenomenon can be explained by the driving effect of large U top and R on the appearance of viscous fingering.
In summary, the increase in U top or R individually is found to suppress density fingering.However, when U top and R grow to a sufficiently high level, they can jointly act as the driving force for viscous fingering.

C. Stability analysis
Following discussions on effects of R and U top , further simulations are undertaken in this section to investigate flow stability.Simulations are conducted for Ra ¼ 10 9 (case III), Ra ¼ 10 8 (case IV), and Ra ¼ 10 7 (case V), with U top ¼ 0-10 À2 and R ¼ À7-7, with simulation parameters being listed in Table II.
To quantitatively illustrate the fingering phenomenon, Fig. 9 demonstrates the l m and t e contour plots in the U top -R plane for cases III when fluid fronts move to l f ¼ 0:75l y .Based on the values of l m , t e , and fingering patterns, the parameter space is divided into five distinct Conversely, in the viscosity-enhanced regime II, a small R enhances density fingering, resulting in large l m and small t e .As U top increases, the displacement effect gradually becomes dominant in the flow, and the parameter space is defined as the displacement-dominated area.
Although U top suppresses density fingering to some extent, it acts as a driving force for viscous fingering.As a result, the gravity-dominated area can be divided into three regimes, namely, viscosity-unstable (III), displacement-suppressed (IV), and stable (V) regimes.In the viscosity-unstable regime III, U top and R are sufficiently large to cause viscous fingering.Consequently, the large l m and small t e are observed in area III.In the displacement-suppressed regime IV, density fingering is suppressed by the moderate U top .In the stable regime V, the large U top inhibits density fingering, while the moderate R is inadequate for inducing viscous fingering.The combination of these two factors results in a stable flow.The fingering patterns corresponding to regions I-V are depicted in Fig. 9(c).
After dividing the parameter space into five regimes according to the flow stability for cases III, similar analysis is conducted to summarize the stable and unstable regions for cases IV (Ra ¼ 10 8 ) and cases V (Ra ¼ 10 7 ).To quantitatively present the flow stability, the fingering length l m =l y at t e (i.e., l f ¼ 0:75l y ) is graphed against R on the x axis, U top on the y axis, and l m =l y on the z axis in Fig. 10  Moreover, a close inspection of the flow regimes in cases III-V reveals the effects of Ra on the flow stability.As Ra decreases, the stable regime V expands in the parameter space, while the gravitydominated regimes I and II contract.This is attributed to the fact that the reduced buoyance force F with the decreasing Ra degrades the density fingering intensity.Specifically, as Ra decreases from 10 9 to 10 8 , density fingering is weakened and thus part of the initially unstable regimes I and II become stabilized.

V. CONCLUSIONS
In this work, a multiple-relaxation-time lattice Boltzmann method is utilized to study the displacements between two miscible and incompressible fluids in porous media at the pore scale.Considering the coexistence of interface instabilities (i.e., density fingering and viscous fingering), simulations are conducted for different values of Rayleigh number Ra, viscosity ratio R, and injected velocity U top .Initially, the general fingering dynamics are discussed in cases without or with U top .From the perspective of fingering development, the flow follows a similar route from the diffusion-dominated stage into the unstable fingering stage.For a fixed Ra, the density fingering intensity is found to be attenuated by the growing R or U top due to the enhanced viscous resistance or displacement magnitude.Consequently, a large U top in conjunction with a moderate R can effectively inhibit the onset of density fingering and introduce a stable flow regime.Notably, with the continuous increase in both U top and R, these two parameters gradually facilitate the development of viscous fingering.In addition, the rise in Ra is shown to amplify the buoyance force and subsequently intensity of the density fingering.These fingering phenomena are investigated quantitatively by evaluating the fingering development time period t e , density jump Dq, front position l f , mixing length l m ; and fingering growth velocity U m .Finally, by comparing values of l m and t e for each fixed Ra, five flow patterns are classified as viscosity-suppressed (I), viscosity-enhanced (II), viscosityunstable (III), displacement-suppressed (IV), and stable (V) regimes.This classification depends on the relative values of U top and R.Moreover, a three-dimensional phase diagram spanned by Ra, U top ; and R is plotted, which demonstrates the distributions and parameter ranges of the above five flow patterns.
Based on the above simulations and discussions, R, U top ; and Ra are proven to have complex impacts on the interface instabilities, thereby influencing the spreading fluids.These findings have profound implications for improving the efficiency of CO 2 capture and storage.Due to the limitation of the 2D system adopted, the porosity of the medium is larger than that of the realistic 3D system.To further advance this field of study, future research should be focused on the mixing process with reactions and realistic porous media in 3D.

SUPPLEMENTARY MATERIAL
See the supplementary material for the validation of the multiplerelaxation-time lattice Boltzmann method for density instability with viscosity contrast.
t ð ÞÀ g j eq x; t ð Þ Â Ã ; (A2) for i; j ¼ 0; 1; …; 8 representing discrete velocity directions, where f i x; t ð Þ and g i x; t ð Þ are distribution functions at location x and time t for the fluid density and the concentration, respectively, e i is the discrete velocity and d t is the time step.The equilibrium distribution functions are defined as 56 where w i is the weights specified to the chosen velocity set, and c s is the lattice sound velocity.q p is a variable related to the fluid pressure as p ¼ c 2 s q p , and q 0 is the fluid density.The distribution function for the force F is defined as The transformation matrix M for the D2Q9 model is This matrix maps distribution functions f and g to the moment space f ¼ M Á f and ĝ ¼ M Á g, respectively.Now, the evolution equations become where S and S C are the diagonal relaxation factor matrices for f and ĝ , respectively.The equilibrium distribution functions in the moment space are described as f eq ¼ q p ; À2q p þ 3q 0 u 2 ; q p À 3q 0 u 2 ; q 0 u; Àq 0 u; q 0 v; Àq 0 v; q 0 u 2 À v 2 ð Þ ; q 0 uv ; (A9) ĝ eq ¼ C 1; À2 þ 3u 2 ; 1 À 3u 2 ; u; Àu; v; Àv; u 2 À v 2 ; uv ð Þ : (A10) The forcing moments F are described as F ¼ 0; 6u Á F; À6u Á F; F x ; ÀF x ; F y ; ÀF y ; 2 uF x À vF y ð Þ ; uF y þ vF x ð Þ : (A11) The macroscopic fluid density, velocity, and concentration can be derived from the distribution functions as Note that here in LB equations, the relaxation times s and s c can be derived from the Chapman-Enskog analysis as 57 As for the boundary conditions, the non-equilibrium extrapolation method is utilized for the top and bottom flows to realize the prescribed velocity and concentration 56 where x f denotes the adjacent location with known distribution functions next to the boundaries.For the no-slip and no-flux boundary conditions at the solid grain surface in porous media, the halfway bounce-back method is utilized as 58 where e i ¼ Àe i , with e i pointing to the solid grains, and f 0 i x; t ð Þ and g 0 i x; t ð Þ representing the post-collision distribution functions.

FIG. 1 .
FIG. 1.The model schematic: two miscible fluids with viscosity contrast and injection velocity in porous media.

FIG. 6 .
FIG. 6.(a) Density contours q Ã when fluid fronts move to l f ¼ 0:75l y , for cases with Ra ¼ 10 9 and different values of U top and R. (b) Plots of the end time t e as a function of R for different U top .

FIG. 7 .
FIG. 7. Plots of (a) the mixing length l m and (b) the growth velocity of fingering U m when l f ¼ 0:75 l y as a function of R for different U top .
. The overall flow instability increases with larger values of l m =l y .The aforementioned five regimes can be identified generally for other values of Ra.In addition, Fig. 10 also confirms the effects of U top and R that the increasing U top or R individually suppresses density fingering, while they together can gradually drive viscous fingering.

FIG. 9 .
FIG. 9. Contours of (a) the mixing length l m and (b) end time t e when fluid fronts move to l f ¼ 0:75l y , in the U top -R plane for case III with Ra ¼ 10 9 .Five distinct regimes are classified as: viscosity-suppressed (I), viscosity-enhanced (II), viscosity-unstable (III), displacement-suppressed (IV), and stable (V).(c) Corresponding density contours are presented to illustrate fingering patterns under five regimes.

TABLE I .
Parameters for case I and II in Part A.