Lattice Boltzmann Simulation of Multicomponent Porous Media Flows With Chemical Reaction

Flows with chemical reactions in porous media are fundamental phenomena encountered in many natural, industrial, and scientific areas. For such flows, most existing studies use continuum assumptions and focus on volume-averaged properties on macroscopic scales. Considering the complex porous structures and fluid–solid interactions in realistic situations, this study develops a sophisticated lattice Boltzmann (LB) model for simulating reactive flows in porous media on the pore scale. In the present model, separate LB equations are built for multicomponent flows and chemical species evolutions, source terms are derived for heat and mass transfer, boundary schemes are formulated for surface reaction, and correction terms are introduced for temperature-dependent density. Thus, the present LB model offers a capability to capture pore-scale information of compressible/incompressible fluid motions, homogeneous reaction between miscible fluids, and heterogeneous reaction at the fluid–solid interface in porous media. Different scenarios of density fingering with homogeneous reaction are investigated, with effects of viscosity contrast being clarified. Furthermore, by introducing thermal flows, the solid coke combustion is modeled in porous media. During coke combustion, fluid viscosity is affected by heat and mass transfer, which results in unstable combustion fronts.


INTRODUCTION
Flows of miscible fluids in porous media are frequently encountered in natural and industrial processes, like underground water movement [1], geologic carbon sequestration [2], enhanced oil recovery [3], and blood transport [4]. Dynamics of such flows are complex due to the unsteady flow, heat and mass transfer, and porous structure involved [5]. This situation becomes more complicated when chemical reactions are introduced.
In a porous medium, two fluids containing chemical solutes are put into contact along an interface. Chemical solutes are transported through the interface to react and thereby modify fluid properties, like density and viscosity. Meanwhile, if density or viscosity gradients exist across the interface, density or viscosity fingering may develop during the fluid transport [6]. Thus, interplay between interface instability and chemical reaction ensues. Such phenomena have been studied, for example, in systems containing ester-alkaline or CO 2 -alkaline miscible solutions. Experiments were carried out to show that changes in chemical solutes modified the reaction type and intensity, leading to either enhanced or suppressed density fingering [7][8][9][10]. In parallel, theoretical analyses were conducted to predict viscous or density fingering scenarios with the reaction A + B → C [11][12][13]. In a porous medium with a partially miscible interface between the two fluids, Loodts et al. [11] suggested that different density contributions and diffusion rates of the three chemical species (A, B, and C) could introduce eight types of density profiles. Trevelyan et al. [12] theoretically showed that such a problem with a miscible interface could yield up to sixty-two types of density profiles. Each profile potentially represented a unique type of density fingering. As for viscous fingering with the reaction A + B → C, a linear stability analysis predicted and classified various fingering scenarios in a parameter plane spanned by viscosity ratios [13]. Under the guidance of theoretical predictions, efforts have been devoted to numerical simulations of interface instabilities with chemical reactions in porous media. For instance, Loodts et al. [14] numerically showed that the reaction A + B → C could accelerate, inhibit, or introduce the development of density fingering. They further demonstrated that different species diffusion rates could introduce four fingering types [15]. Numerical simulations were also performed to model viscous fingering dynamics with the reaction A + B → C [16]. Results suggested that the reaction could introduce a non-monotonic viscosity profile and thereby trigger the development of fingering. These simulations provided spatial and temporal properties of interface instabilities with chemical reactions, thus enriching experimental and theoretical findings.
The above studies focused on porous media flows with a homogeneous reaction between miscible fluids. On the other hand, a heterogeneous reaction can take place between fluid and solid phases. Explicitly, such a type of reaction consumes chemical solutes from the fluid phase and solid matrices to yield dissolved (dissolution reaction) or solid (precipitation reaction) products, which subsequently modifies medium structures and flow motions [17]. Until now, porous media flows with heterogeneous reactions have been extensively investigated. For example, acid dissolution experiments showed that the dissolution reaction increased medium porosity behind the moving fluid interface and thereby triggered fingering phenomena (or infiltration instability) [6,18]. Fredd and Fogler were among the first to experimentally identify different fingering channels under various flow and dissolution rates [19]. The fingering phenomena induced by the dissolution reaction were observed to be destabilized by medium heterogeneity and brine pH [20]. In parallel, experiments on porous media flows with precipitation reactions showed that the reaction could generate a locally adverse mobility gradient to cause the onset of fingering [21]. The interplay between viscous fingering and precipitation reaction was experimentally studied, with the influencing factors and fingering shapes being clarified [22]. In the meantime, numerical simulations of porous media flows with heterogeneous reactions were conducted. For instance, fingering dynamics induced by dissolution reactions were tracked and compared on different simulation length scales [23]. Continuum models were built to describe fingering properties and classify different fingering regimes under the effects of chemical dissolution [24,25]. In the presence of precipitation reactions, viscous fingering was modeled and found to be intensified [26,27]. In some applications, like in situ combustion for heavy oil recovery, the combustion between hot air and solid fuels, as a typical form of heterogeneous reaction, releases heat to change fluid properties [28]. Thus, thermal flows were included in numerical studies on porous media flows with heterogeneous reaction. Field-scale modeling of solid coke combustion suggested the important role of heat conduction in stabilizing the combustion front [29] and identified key factors influencing the coke combustion rate and the front sweep efficiency [30,31].
Overall, for porous media flows with chemical reactions, existing works focus on using empirical correlations and volume-averaged techniques to determine effective flow and reaction properties on microscopic scales. However, pore-scale information is necessary to provide further understanding of the microscopic phenomena involved [32]. In the past 3 decades, the lattice Boltzmann (LB) method has been developed and has become a powerful solver for porous media flows on the pore scale. This is attributed to its attractive advantages, like its general kinetic foundation, high parallelism, and ability to handle multiphysics and complex boundaries [33,34]. Recently, LB models have been proposed to simulate flows of reactive fluids in porous media, like methane hydrate dissolution [35], solid coke combustion [32], and reactive mixing with viscous fingering [36]. Although these works showed the superior ability of the LB method in modeling porous media flows, some limitations should be noted. For example, changes in fluid density or viscosity with temperature are usually not accounted for, models for multicomponent reactions (like A + B → C) are underdeveloped, and iterative boundary schemes for conjugate heat transfer are complex and time-consuming. We have developed a series of LB models to overcome these shortcomings separately [37][38][39][40]. Nevertheless, a uniform LB model is still needed. This work thus proposes an LB model to simulate pore-scale reactive flows in porous media, with homogeneous and heterogeneous reactions and compressible and incompressible fluid densities being considered.

MATHEMATICAL MODEL
In a porous medium with porosity ϕ, a fluid 2 of solute B initially fills the pore spaces between solid matrices. Another fluid 1 of solute A is put into contact with fluid 2 along an interface. The two fluids are assumed to be miscible and have their own species concentration, viscosity, temperature, and density, which are noted by subscripts 1 and 2, respectively. As an example, a schematic of the fluid distribution with a vertical interface is shown in Figure 1. Once fluid transport starts, the two solutes A and B come into contact and react in the mixing zone as follows: The reaction product C is in the dissolved form, Q is the reaction heat, and the reaction rate of this homogeneous reaction is R o . Meanwhile, the nonreactive solid matrices can be coated with a reactive solid Bs to react with A as follows: The product P is in either the dissolved or the solid form, representing the dissolution or the precipitation reaction, respectively. It should be noted that fluid 2 does not contain solute A, and thus, this heterogeneous reaction only takes place at the interface I between fluid 1 and solid Bs. With this heterogeneous reaction, solid Bs is dissolved gradually, and the update of the porous structure is tracked as follows [41]: where V B and V Bm are the volume and the molar volume of solid Bs, respectively, and S B is the reactive surface area. The reaction rate of this heterogeneous reaction is noted as R e . The above homogeneous and heterogeneous reactions usually cause variations in fluid properties (like concentration, temperature, density, and viscosity) and medium structures, which subsequently alter flow motions in porous media. We seek to study the fluid transport properties of such a system, focusing on effects of chemical reactions. Evolutions of fluid flow and species concentrations in pore spaces are described by the following conservation equations: where u (u, v), ρ l , p, and ] are the fluid velocity, density, pressure, and kinematic viscosity, respectively. t is the time, and F is the body force. Y r is the mass fraction of species r (r A, B, C, P), and it is related to species concentration as C r ρ l Y r /M r , with M r being the species molecular weight. D r and S r are the diffusion coefficient and the homogeneous reaction source term of species r (Eq. 1), respectively. Furthermore, chemical reactions in Eqs 1, 2 usually release heat to change the local temperature, which in turn modifies the reaction rate and fluid properties. For such a consideration, heat transfer in both the fluid and solid phases is built as follows: Here T, c p , α, and ρ are the local temperature, specific heat at constant pressure, thermal diffusivity, and density, respectively. In addition, at the interface I between solute A and solid Bs, the heterogeneous reaction Eq. 2 is described as follows [32,41]: Species conservation : Conjugate heat transfer : T I,+ T I,− , n · k∇T + ρc p uT where n is the interface normal pointing to the fluid phase, + and − denote parameters on either side of I, k αρc p is the local thermal conductivity, q is the reaction heat flux, and a r is the stoichiometric coefficient of species r. By introducing the characteristic length L, velocity U, temperature T ch , concentration C ch , and density ρ ch , dimensionless parameters marked by asterisks are derived as follows: Key characteristic numbers are as follows: the Reynolds number Re, the Peclet number Pe, the Prandtl number Pr, and the Schmidt number Sc.
It is emphasized that the above equations consider homogeneous and heterogeneous reactions via the source term S r and boundary conditions Eqs. 8-10, respectively. Evolutions of chemical species are described using Eq. 6 separately, thus allowing for homogeneous reactions between multiple solutes. Furthermore, both compressible and incompressible fluid flows can be taken into account, by setting the fluid density as a constant and a variable, respectively. On one hand, for compressible flows, the inclusion of heat transfer brings in the temperature-dependent fluid density. On the other hand, for incompressible flows, the thermophysical properties in Eq. 7 are fixed as constants in the fluid phase ϑ l ϑ 0 (ϑ c p , α, ρ), which can be different from those in the solid phase, namely, ϑ s ≠ ϑ l .

LATTICE BOLTZMANN MODEL
An LB model is now developed to solve the above governing equations for porous media flows with chemical reactions. For pore-scale simulations, the multiple-relaxation-time (MRT) scheme was shown to be able to correctly capture the no-slip velocity condition at the fluid-solid boundary and avoid the unphysical dependence of permeability on viscosity [42][43][44]. The MRT collision operator is thus employed in our model development.
Considering that this work focuses on reactive fluid flows in twodimensional (2D) porous media, the most popular two-dimensional nine-velocity (D2Q9) scheme is applied, with the discrete velocities e i and weight coefficients w i being as follows [34]: Here, e δ x /δ t is the lattice speed, with δ x and δ t denoting the lattice spacing and the time step, respectively. In this work, e is given as the velocity unit, that is, e 1. The transformation matrix M corresponding to the D2Q9 scheme is as follows [34]: This matrix maps the distribution functions from the physical space ψ (ψ 0 , ψ 1 , ψ 2 , . . . , ψ 8 ) T into the moment space as follows [34]: Due to the temperature-dependent density and different thermophysical properties between the solid and fluid phases, Eqs. 6, 7 are rewritten as follows [45]: Source terms in the derived species and energy conservation equations are as follows [39]: More details about the derivation can be found in our recent work [39]. It should be noted that that for compressible fluids with fixed ρ l ρ 0 , the terms F r2 and z t (ρc p ) in F T2 are reduced to zero. An MRT LB model is developed to solve Eqs 4, 5, 15, 16 for describing reactive fluid flows in porous media. The LB evolution equations are as follows [34,37,46]: , and h i (x, t) are distribution functions for the fluid density, mass fraction of species r, and local temperature, respectively. The corresponding equilibrium distribution functions f eq i , g eq r,i , and h eq i are given as follows [34,46,47]: Here, ρ p is a variable related to the fluid pressure as p c 2 s ρ p , with c s e/ 3 √ being the lattice sound velocity. Different expressions of f eq i are applied to account for compressible and incompressible fluid flows, respectively. The correlation term Λ is given as follows: with D and θ being the model dimension and the dimensionless temperature of a fluid mixture, respectively. The other correction term C̄j in Eq. 18 is introduced to eliminate the deviation from third-order velocity moments [46]. With these two modifications, thermal expansion effects in compressible fluid flows are realized via the temperature-dependent density. The equation of state to connect the fluid pressure and the temperature is defined as p ρrT. On the other hand, for incompressible fluid flows, these two correlation terms are reduced to zero, and the dimensionless temperature becomes θ 1. In evolution Eqs. 18-20, S, S r , and S T are diagonal relaxation matrices of the relaxation rates s i , s r,i , and s T,i in the moment space, respectively. To avoid discrete lattice effects, distribution functions for the force and source terms are defined as follows [34,36]: Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 715791 Time derivative terms in Eqs 19, 20 are treated using the backward-difference scheme as follows [5]: Based on the transformation equation (14), evolution equations (18)(19)(20) are implemented in the moment space as follows: ĥx The equilibrium moments fêq, ĝe q r , and hêq are expressed as follows: and the corresponding moments for the force and source terms are as follows: F̂r F r 1, −2, 1, 1 − 0.5s r,3 u, − 1 − 0.5s r,4 u, Finally, macroscopic variables are obtained from the distribution functions as follows: Through the Chapman-Enskog analysis on the present LB equations, the governing equations can be recovered with the relaxation times τ, τ r , and τ T being as follows [46,47]: Furthermore, for compressible fluid flows, the constraints on the correction term in the moment space (Ĉ M · C) are formulated as follows [46]: with Q x ρu(1 − θ − u 2 ) and Q y ρv(1 − θ − v 2 ). By such an analysis, the gradient terms of the mass fraction and temperature in Eq. 17 are determined as follows [36]: The other gradient terms in Eq. 17 are calculated using the isotropic central scheme as follows [48]: Finally, the heterogeneous reaction at interface I between fluid 1 and solid Bs is modeled by resolving the interface conditions in Eqs. 8-10. On one hand, with the derived thermal source term F T , conjugate heat transfer conditions (Eq. 10) are realized by solving the energy Eq. 7. On the other hand, to solve Eqs 8, 9, the interface mass fraction gradient ∇Y I r is first determined using the finite-difference scheme as follows [49]: where Y l r is the species mass fraction at the fluid node neighboring I. By inserting Eq. 43 into Eq. 9, the value of Y I r is determined. Thus, Eqs. 8, 9 describe a reactive boundary with no-slip velocity and a given Y I r , which are implemented via the halfway bounce-back scheme. The unknown distribution functions at the fluid node x f adjacent to I are as follows [49]: g r,ı̄x where e ı̄ −e i , with e i pointing to the solid phase. f i ′ (x f , t) and g r,i ′ (x f , t) denote post-collision values.
Compared with recent LB models for porous media flows [32,35], the present MRT LB model offers several advances. First, by introducing correction terms, both incompressible fluid flows with a fixed fluid density and compressible flows with a temperaturedependent one can be modeled. Second, separate LB equations with reactive source terms are developed to account for homogeneous reactions between multiple solutes, like A + B → C. Third, for heterogeneous reaction at the solid-fluid interface, additional source terms and the bounce-back scheme are built to realize the conjugate heat transfer and species conservation conditions, without iterative calculations. In addition, different from our recent LB models [37][38][39][40], the present one offers a capacity to achieve these advances simultaneously.
In the above section, a multicomponent MRT LB model has been presented for studying porous media flows with chemical reactions on the pore scale. Advances of this model over existing ones are introduced. This section conducts simulations of reactive fluid flows in porous media. The obtained results will show the ability of the present LB model in predicting hydrodynamic instabilities at the fluid interface, thermal flows in both pores and solid matrices, and effects of chemical reactions on fluid transport.

Density Fingering With the Reaction
This section aims to investigate density fingering with the reaction A+ B → C between two miscible solutions 1 and 2 in porous media. As depicted in Figure 2, a 2D homogeneous structure is constructed, with the geometric parameters l x 1, l y 2/3, r x 27δ x , r y 27δ x , d 12δ x , l p 15δ x , and ϕ 0.69 (in lattice units). Initially, the porous medium is saturated with a host fluid 2 of solute B. Another fluid 1 of solute A is placed upon the medium and in contact with fluid 2 along the top boundary (y 0). The two solutions are assumed to be incompressible and isothermal. Then, across the partially miscible top boundary, A dissolves into fluid 2, and no mass transfers in the reverse direction. Dissolved A reacts with B to generate a solute C as in Eq. 1, with the reaction rate being as follows [14]: Here, k is the kinetic reaction constant. All three chemical species contribute to fluid density and viscosity changes as follows [13,14]: where β r is the concentration expansion coefficient of species r. R A ln(μ A /μ B ) and R C ln(μ C /μ B ) are log-to-viscosity ratios between pure solutions of species r at concentration C r,0 . Under the gravity field, the body force (or buoyancy force) is calculated as follows [37]: Here, g is the acceleration vector of gravity. With such effects, different density stratifications develop, and density fingering may be triggered. In this incompressible system, the thermal flow is not considered. Thus, fluid motion and species evolutions in pore spaces are described using Eqs 4-6. In addition, fluid density is fixed as ρ l ρ 0 , the reaction source term is defined as S r a r R o , the heterogeneous reaction at the fluid-solid interface is not considered, and the three species are assumed to diffuse at the same rate D r D. For dimensionless parameters in Eq. 11, the characteristic length L, velocity U, concentration C ch , and density ρ ch are as follows: The Rayleigh numbers Ra r and the Damköhler number Da are defined as follows: For such a problem, our recent works have studied fingering dynamics with effects of the reaction A + B → C and differential diffusion [37,38]. Thus, in the following simulations, we introduce the impact of the three chemical species on fluid viscosity (Eq. 48) and seek to understand how the development of density fingering is affected. The Schmidt number and the Damköhler number are fixed as Sc 100 and Da 5, respectively, and different values of Ra AB , Ra CB , R A , and R C are selected to change test conditions. During LB simulations, the no-slip and no-flux bottom wall and matrix interface are realized using the halfway bounce-back scheme, periodic conditions are set at the two lateral boundaries, and the partially miscible top boundary is implemented using the nonequilibrium extrapolation scheme. The relaxation rates in this simulation case are set as follows: s 0 s 3 s 5 0, s 1 1.1, s 2 1.2, s 4 s 6 (16τ − 8)/(8τ − 1), s 7 s 8 1/τ, s r,0 1, s r,1 s r,2 1.1, s r, 3 s r,4 s r,5 s r,6 1/τ r , and s r,7 s r,8 0.9. A mesh of N x × N y 1,500 × 1,000 lattice size is used after grid convergence validations.
We focus on density fingering phenomena with fluid viscosity variations caused by chemical species. Simulations with Ra AB 1, Ra CB 1, R A 0, and different values of R C are first performed. From our calculated results, different types of density fingering are identified, and fingering is found to be stabilized as C decreases the fluid viscosity (R C < 0). To illustrate such effects, the results of three cases with R C −3, 0, 3 are provided and discussed. It should be noted that the three values of R C represent the three different scenarios where the reaction product C decreases, cannot change, or increases the fluid viscosity, respectively. For each of the three cases, Figure 3 depicts density distributions at six time instants, with the dimensionless density being calculated as ρ* (ρ − ρ 0 )/ (ρ 0 β B C B,0 ). From the constant-viscosity case (R C 0) in Figures 3A,B, a classical density fingering development pattern is observed, that is, fluid diffusion, individual fingering growth, fingering merge, and new fingering reinitiation [37]. It should be noted that during the initial diffusion stage, dense fluid accumulates near the top boundary to overcome the viscous shear force blow. Fingering develops when the driving force introduced by the accumulated dense fluid is large enough.
Compared with this constant-viscosity case, the introduction of viscosity variations in Figures 3A,C yields different density fingering dynamics. On one hand, when chemical product C decreases the fluid viscosity ( Figure 3A with R C −3), the classical fingering development pattern is observed. However, fingering starts earlier and grows faster than in the constant-viscosity case ( Figure 3B). Such an accelerated fingering propagation finally leads to diluted fingers. On the other hand, in the case with solute C increasing the fluid viscosity ( Figure 3C with R C 3), even though the dissolution of A and the reaction A+ B → C can introduce a buoyantly unstable stratification, the fluid interface remains flat at each time instant, and no density fingering occurs.
In order to explain such influence of fluid viscosity variations, distributions of species concentrations (C * r C r /C B,0 ) and fluid viscosity (μ* μ/μ 2 ) are calculated and provided in Figure 4. For each case, results in Figure 4 correspond to the density field in Figure 3 at the last time instant. First, in the constant-viscosity case, Figure 4B shows that μ* distributes uniformly and density fingering is introduced by fingers of species A. Then, in the case with R C −3 ( Figure 4A), fingers of A are also observed to be at the origin of density fingering. However, compared with the constant-viscosity case ( Figure 4B), μ* (or fluid viscous force) decreases with the generation of C in the top layer. It implies that less A is required to accumulate near the top boundary to overcome the viscous force. Thus, fingers of A start earlier and develop faster in the top layer occupied with C. Finally, in Figure 4C with R C 3, the chemical product C increases the fluid viscosity and viscous force in the top area. Therefore, even  though dissolved A accumulates to overcome the enlarged viscous force, the driving force introduced by the accumulation of A is still not large enough, and the fluid interface remains flat without fingering phenomena. These results explain the observations in Figure 3.
To quantitatively compare fingering dynamics in the three cases, Figure 5 plots temporal evolutions of the front position l m /l y , volume-averaged storage of A in the host fluid 〈C * A + C * C 〉, and horizontally averaged mass flux of A at the top boundary J*. Here, l m is defined as the most advanced vertical position of density fingering in the host fluid, and J* is calculated as ). Curves of l m /l y show that the fingering front keeps growing in a diffusive tendency at R C 3, indicating the stable fluid interface. In the other two cases, however, l m /l y gradually departs from the diffusive tendency due to the onset of density fingering. Moreover, among the three cases, fingering at R C −3 propagates the fastest. This is because solute C decreases the fluid viscosity and destabilizes fingering development. Results of 〈C * A + C * C 〉 and J* also suggest that the decreased fluid viscosity at R C −3 accelerates the storage speed of A in the host fluid. Finally, by comparing the status of the three cases with l m /l y 0.5, the stored amount of A is found to be the smallest in the case R C −3. This is attributed to the fact that the fast fingering propagation decreases the time period for the dissolution of A from the top boundary. These results quantitatively confirm the above fingering scenarios in Figure 3.
In general, when solute C decreases the fluid viscosity (R C < 0), density fingering is destabilized, and the destabilizing intensity ascends with descending R C . To further illustrate this impact, three cases with Ra AB 0.1 are conducted. Different from the above cases with Ra AB 1, the contribution to fluid density of A is smaller than that of B. Thus, in nonreactive cases, the unstable density contrast introduced by the dissolution of A from the top boundary is not large enough to bring in density fingering phenomena. As conducted in cases with Ra AB 1, Figures 6,  7 provide distributions of fluid density, species concentrations,  and fluid viscosity. In the constant-viscosity case with Ra CB 1 and R C 0 ( Figures 6B, 7B), the uniform fluid viscosity distribution and the flat fluid interface are observed. This indicates that the reaction A + B → C, without viscosity variations, cannot trigger the onset of density fingering. On the contrary, once the unstable viscosity contrast (R C −3) comes into play in the other two cases, density fingering appears and experiences the classical fingering development stages. First, in the case with Ra CB 1 and R C −3 ( Figures  6A, 7A), the density contrast remains the same as in the constantviscosity case ( Figures 6B, 7B), but fluid viscosity in the top layer drops down with the generation of C. As a result, viscous force decreases and fluid mobility increases in the top area, which makes dissolved A grow in a finger-like shape and thus induces the development of density fingering. Second, in the case with Ra CB 0.1 and R C −3 ( Figures 6C, 7C), a buoyantly stable density stratification develops with chemical reaction. Nevertheless, as illustrated in the zoom-in view, the fluid interface is distorted with time and finally becomes fingered. It is because the chemical product C decreases fluid viscosity in the top area and thus helps dissolved A develop into density fingering. These cases with a small density contrast suggest that decreased fluid viscosity with the chemical product C can destabilize or even trigger the development of density fingering.

Solid Coke Combustion With Viscous Fingering
The developed LB model is then applied to explore the solid coke combustion dynamics in porous media, which is a typical heterogeneous reaction at the interface between hot air and coke. As displayed in Figure 8, a 2D homogeneous structure with porosity ϕ 0.406 is considered. The geometric parameters are set as l x 900 μm, l y 760 μm, d 18 μm, δ 3 μm, r x 42 μm, and r y 33 μm, respectively. Initially, solid coke is homogeneously deposited on unreactive solid matrices, and the pore spaces are filled with hot nitrogen at the temperature T 0 773 K. Then, hot air at the temperature T 0 and with the oxygen (O 2 ) mass fraction Y O2,0 0.233 is injected to react with coke by a driving force F (F x , 0). For simplicity, coke combustion is assumed to take place at the air-coke interface I as C + O 2 → CO 2 + Q [50]. The released heat is calculated using Q R e h r , with h r being the chemical reaction heat. The reaction rate R e is estimated according to the first-order Arrheniusttype Eq. 39 as follows: where A, E, and R are the pre-exponential factor, the activation energy, and the ideal gas constant, respectively. As for boundary conditions, zero-gradient conditions are applied for all the scalars at the outlet, and the top and the bottom are periodic. The impact of reaction heat and mass transfer on fluid viscosity is considered as follows [51]: Here, R y ln(μ 1 /μ 2 ) and R t ln(μ T h /μ 2 ) are the mass and thermal viscosity ratios, respectively. μ 1 and μ 2 are viscosities of fluids 1 and 2 at T T 0 , respectively. μ T h is the viscosity of fluid 2 at T T h .
For this case, fluid flow and species evolutions in pore spaces, as well as heat transfer in both pores and solid phases are described using Eqs 4-7. The reaction source term is set as S r 0 because homogeneous reaction is not involved. On the other hand, coke combustion is described using Eqs 8-10 and  implemented using the bounce-back scheme in Eqs 44, 45. The conversion between physical and lattice units is based on dimensionless parameters in Eq. 11 and the characteristic parameters are selected as L 0.5r x , U α l /L, ρ ch ρ l,0 , and T ch 3,000 K [39]. In the subsequent simulations, the required parameters are set as A 9.717 × 10 6 m/s, E 131.09 kJ/mol, hr 388.5 kJ/mol, ρ s c p,s /ρ l c p,l 368, and α s /α g 0.119, respectively. The Reynolds, Prandtl, and Peclet numbers are fixed as Re 1.4, Pr 0.7, and Pe 1.5, the driving force is F * x 7.3, and a mesh of size N x × N y 900 × 760 is used after grid convergence tests. In addition, the relaxation rates used in this case are as follows: s 0 s 3 s 5 1, s 1 s 2 0.9, s 4 s 6 (16τ − 8)/(8τ − 1), s 7 s 8 1/τ, s r,0 1, s r,1 s r,2 1.1, s r,3 s r,4 s r,5 s r, 6 1/τ r , s r,7 s r,8 0.9, s T,0 1, s T,1 s T,2 1.1, s T,3 s T,4 s T,5 s T,6 1/τ T , and s T,7 s T, Our recent work has investigated the general coke combustion dynamics and evaluated the key influencing conditions, with the fluid viscosity being assumed to be a constant [39]. In practical applications, however, fluid viscosity varies with heat and mass transfer. Therefore, this work extends to studying the impact of viscosity variations on the combustion front stability. A base case with R y 0 and R t 0 is simulated at first. Figure 9 illustrates the distributions of residual coke, O 2 mass fraction Y O 2 , and temperature T/T 0 at the time instant t 5.62 s. As can be seen, the combustion front extends around a grain diameter. In this area, coke and O 2 are fully consumed by combustion, and the combustion front remains flat without fingering phenomena. Furthermore, temperature is observed to increase from the inlet and then reaches a peak value at the combustion front and finally drops to T 0 on the downstream side. This distribution is attributed to the released combustion heat at the combustion front. It is also noticed that fluid on the downstream side (i.e., the area from the combustion front to the outside) is slightly hotter than that on the upstream side. This is caused by both the flow direction and the heat transfer that is faster than the combustion front movement [39]. In this case, the coke combustion with the stable combustion front, low temperature increase, and full coke and O 2 utilization rate can fit engineering requirements.
The fluid viscosity remains uniform in the above base case. Thus, two cases with R y −4, R t 0 and R y −4, R t −3 are then considered. Other simulation parameters are set to be the same as in the base case, except for the Peclet number Pe 10. Spatial distributions of residual coke, the O 2 mass fraction Y O 2 , and temperature T/T 0 for the two cases are provided in Figure 10. Results show that in each of the two cases with R y −4, the combustion front proceeds in a finger-like pattern along the x direction. This is different from the flat combustion front in the base case (Figure 9), implying the destabilizing effects of the viscosity contrast R y −4. As for the temperature field, the maximum temperature occurs at the combustion front as in the base case. However, compared with the base case, the combustion fronts in these two cases are decelerated, and the temperature of the downstream side is no longer obviously higher than that of the upstream side. This is because the highly viscous fluid on the downstream side prevents the fluid propagation and thus slows down both the front movement and the heat transfer.
In the case with R t −3, the fluid viscosity decreases with the increasing temperature. By comparing results in Figures 10A,B, such temperature-induced viscosity variations are found to influence coke combustion properties. On one hand, the fluid prorogation along the x direction is accelerated and more fingers at the combustion front are triggered. This destabilizing impact is attributed to the fact that combustion heat is released to increase the local temperature and decrease the fluid viscosity at the combustion front. Subsequently, the unstable viscosity contrast between fluids at the combustion front and on the downstream side is enlarged, which enhances viscous fingering phenomena. On the other hand, the combustion front temperature increases. It is explained by the fact that in the combustion front area, the high temperature induces the low fluid viscosity and the accelerated fluid flow velocity. Subsequently, the injected O 2 can FIGURE 9 | Contours of residual coke, the O 2 mass fraction Y O2 , and temperature T/T 0 for the base case with R y 0 and R t 0 at t 5.62 s. FIGURE 10 | Contours of residual coke, the O 2 mass fraction Y O2 , and temperature T/T 0 for the two cases at t 8.14 s with (A) R y 4, R t 0 and (B) R y 4, R t −3.
Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 715791 easily expand to react with coke and release heat, which leads to the increased combustion front temperature.
To further quantify differences between the two cases with viscosity variations, temporal evolutions of the volume-averaged temperature T v /T 0 and the residual coke ratio Vr c are calculated and illustrated in Figure 11. As can be seen, each line of Vr c decreases with time because coke burns out gradually, while every temperature curve shows an increasing pattern due to the released combustion heat. In addition, the inclusion of temperature-induced viscosity variations is found to bring in high burning temperature and fast coke consumption. It is due to the fact that the increased front temperature brings in the decreased fluid viscosity and, hence, the enhanced coke combustion intensity. This quantitatively verifies the above observations in Figure 10. It should be emphasized that these two cases feature finger-like combustion fronts. They are thus undesirable in practical applications because the downstream fluid cannot be efficiently displaced.

CONCLUSION
Flows of reactive fluids in porous media have been studied on the pore scale. A multicomponent multiple-relaxation-time lattice Boltzmann (LB) model is developed to describe fluid motions and species evolutions in pore spaces, as well as heat transfer in both void pores and the solid phase. In the present model, separate LB equations are built to describe fluids and species evolutions during homogeneous reaction between miscible fluids; source terms and boundary schemes are derived to account for heterogeneous reaction at the fluid-solid interface; and additional correlation terms are introduced to model both incompressible and compressible (with temperature-dependent fluid density) fluid flows. Based on this model, two types of physical problems are simulated. One is density fingering with the homogeneous reaction A + B → C. Results show that depending on variations in the fluid viscosity caused by the chemical product C, the development of density fingering can be enhanced, suppressed, or even triggered. To further study effects of thermal flows, solid coke combustion in porous media is carried out. Simulation results successfully capture coke combustion dynamics and the evolution of the porous structure. Furthermore, the viscosity contrast caused by heat and mass transfer is found to change the coke combustion stability considerably. These two cases verify the capability of the proposed LB model in resolving flows of reactive fluids and the coupled viscosity and density instabilities on the pore scale. Extension to 3D porous media flows with chemical reaction is currently underway. The results are expected to clarify the differences and similarities between 2D and 3D simulations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.