Developing energy storage systems or technologies can provide long-term support for future low-carbon energy systems while reducing energy supply risk. At present, the main energy storage technology is pumped hydro energy storage (Rehman et al., 2015; Javed et al., 2020), and the research and application of phase change energy storage (Yang et al., 2020; Huo et al., 2022; Yang et al., 2022; Liu et al., 2023a), battery energy storage (Heyhat et al., 2020; Naghavi Sanjani et al., 2023) and other energy storage technologies are developing (Koohi-Fayegh and Rosen, 2020; Liu et al., 2023c). However, the energy storage technologies mentioned above are still unable to meet the demands for big capacity and long-term energy storage. Meanwhile, underground energy storage (Strobel et al., 2020; Zivar et al., 2021) can serve as both energy transmission and storage in the energy market, and it is a viable solution to the problem of big-capacity long-term energy storage.

Hydrogen energy (Lackey et al., 2023; Zhang et al., 2023) is an excellent solution to the problem of energy sustainability, with advantages such as a large number of sources, a high calorific value, no pollution, and a wide range of applications. Hydrogen has a lower molecular weight than natural gas, hence it requires more space and better sealing. Reservoir rocks, as a result, provide a conducive environment for large-scale hydrogen energy storage and are potential hydrogen storage sites (Lankof and Tarkowski, 2020). The process of underground hydrogen storage is extremely complicated because it involves the interaction of fluid flow with heat and mass transfer in porous media. One common issue in porous media is miscible displacement.

Miscible displacement refers to the displacement of two or more miscible fluids, such as seawater and freshwater, surface sewage and groundwater, tracer-containing fluids and pure fluids, supercritical carbon dioxide and crude oil, which are widely used in the fields of hydrology (Tosco et al., 2014), chemistry, medicine, and petroleum engineering (Jia et al., 2019; Bashir et al., 2022). Due to the differences in the pore structure characteristics (connectivity, tortuosity, etc.) and fluid flow characteristics (viscosity ratio of displacing fluid and displaced fluid, molecular diffusion, displacement flow rate, etc.) of porous media, fingering instabilities often occur during miscible displacement. When a less viscous fluid is intruded into a more viscous fluid, the interface between them can become unstable, resulting in viscous fingering (Saffman and Taylor, 1958). The main focus of this investigation is on the viscous fingering in the miscible displacement process.

The viscous fingering phenomenon has been extensively studied over the past several decades (Homsy, 1987; Bacri et al., 1991; Liu et al., 2023b), with the majority of the study being done at the representative elementary volume (REV) scale (Zimmerman and Homsy, 1992; De Wit and Homsy, 1997; Norouzi and Shoghi, 2014). At the REV scale, a porous medium containing a complex solid skeleton is averaged to a homogeneous medium. The structural properties of the porous medium are characterized only by two macroscopic structural parameters of the porous medium, namely, porosity and permeability. Comparatively, there is relatively less work focused on pore-scale viscous fingering. Viscous fingering involves complex processes such as fluid flow, diffusive mass transfer, and interfacial instabilities. Furthermore, at the pore scale, the displacement process in different directions of a porous medium significantly influences the development of viscous fingering and the progression of the displacing fluid. These processes are closely coupled with pore structure, and the development of the viscous fingering between large and small pores greatly determines the displacement process. Thus, the pore-scale investigation is important and can largely improve our understanding of the effects of these coupling processes and the microscopic pore structure on viscous fingering and displacement efficiency.

Existing studies on the miscible viscous fingering in porous media at the pore scale have mostly focused on the influence of factors such as viscosity ratio, porous media structure (Liu and Guo, 2015; Elgahawy and Azaiez, 2021), and chemical reactions (Lei and Luo, 2019; 2021), ignoring the effect of the gravity field. However, the interplay between the viscosity ratio and the gravity determines the interfacial instabilities (Jiao and Hötzl, 2004; Jiao and Maxworthy, 2008; Zeeshan Mohiuddin and Stokes, 2013). The viscosity ratio and the gravity either stabilize or destabilize the interface. The instabilities of the miscible displacement interface in porous media are determined by the angle of fluid flow relative to the direction of gravity and the viscosity ratio.

In this study, we use the lattice Boltzmann method to simulate viscous fingering in an inclined porous medium under miscible conditions. The complicated coupling between fluid flow, mass transfer, and pore structure is investigated. The effects of the inclination angle, viscosity ratio, and the structure of the porous media on displacement efficiency are also discussed.

In this study, the grayscale image of the porous medium is obtained by using computed tomography (CT) scanning technology, and then the two components of the porous medium skeleton and porous medium pores in the grayscale image are represented by different gray values, to obtain a digital porous medium. The definition of the porosity of the porous medium is φ = V pore V total , (1) where V _pore represents the volume occupied by the pores, and V _total is the total volume of the porous medium.

We then simulate the displacement between two fluids in the two-dimensional porous medium of length L and width W, with gravity acting in the vertical direction, as shown in Figure 1. The porous medium is inclined at an angle θ to the horizontal, the x and y-axes are taken along and perpendicular to the porous medium, respectively. The porosity of the porous medium is φ = 0.865. At first, the displaced fluid (fluid 2) with the kinematic viscosity ν ₂ occupies the porous medium, and then the displacing fluid (fluid 1) with the kinematic viscosity ν ₁ is injected from the left boundary, the average fluid velocity along the x direction is maintained constant and equal to u _in . The upper and lower boundaries of the system are no-slip and no-flux boundaries. Meanwhile, the solid skeleton is impermeable. Therefore, no-slip and no-flux boundaries are also adopted at the fluid-solid interface. At the outlet boundary, the Neumann boundary condition is used.

Assuming that the two fluids are incompressible, the Boussinesq approximation can be adopted. The mixture density ρ and viscosity μ are assumed to be ρ = ρ 2 1 − β C C , (2) μ C = μ 2 e − R C , (3) where C ∈ [0, 1] is the mixture concentration, which represents the volume fraction of the displacing fluid; ρ ₂ is the density of the displaced fluid, β _C is the solute expansion coefficient; μ ₂ is the viscosity of the displaced fluid, R = ln M is the natural logarithm of the viscosity ratio (Tan and Homsy, 1986; Kuang et al., 2003), M = μ ₂/μ ₁ is the viscosity ratio of the displacing and displaced fluid. It is assumed that the fluids’ diffusivity, D, is constant and independent of concentration.

The governing equations in this study include the continuity and Navier-Stokes equations, and the convection-diffusion equation for the concentration of the displacing fluid. The dimensionless governing equations can be written as follows: ∂ U ∂ X + ∂ V ∂ Y = 0 , (4) ∂ U ∂ t + U ∂ U ∂ X + V ∂ U ∂ Y = − ∂ P ∂ X + ∂ ∂ X ν C ∂ U ∂ X + ∂ ∂ Y ν C ∂ U ∂ Y + R a S c C − C 0 sin ⁡ θ , (5) ∂ V ∂ t + U ∂ V ∂ X + V ∂ V ∂ Y = − ∂ P ∂ Y + ∂ ∂ X ν C ∂ V ∂ X + ∂ ∂ Y ν C ∂ V ∂ Y + R a S c C − C 0 cos ⁡ θ , (6) ∂ C ∂ t + U ∂ C ∂ X + V ∂ C ∂ Y = 1 S c ∂ 2 C ∂ X 2 + ∂ 2 C ∂ Y 2 . (7) The dimensionless parameters in Eqs 4–7 are defined as follows, X = x x * , Y = y y * , x * = y * = L * , (8) U = u u * , V = v u * , t = τ t * , P = p p * , (9) R a = g β C Δ C L * 3 ν 2 D , S c = ν 2 D , (10) where u, v represent the velocity component in x and y directions, respectively; L* is the characteristic length, where L* is the width of the porous medium W; τ and p denote time and pressure; u* = ν ₂/W, t* = W ²/ν ₂, p * = ρ ( ν 2 / W ) 2 ; X and Y are dimensionless horizontal and vertical coordinates, respectively; U and V are dimensionless velocity components along X and Y directions, respectively; t and P are the dimensionless time and pressure; ν ₂ is the kinematic viscosity of the displaced fluid; ν(C) = e ^−RC is the dimensionless kinematic viscosity of mixture; Ra is the concentration Rayleigh number, Sc is the Schmidt number; C ₀ = 0.5, g is the gravitational acceleration, and ΔC = 1 is the concentration difference.

In this study, the lattice Boltzmann Method (LBM) is employed to solve the governing equations. Two distribution functions are used to simulate the velocity field and the concentration field, respectively.

The velocity field is described by the following evolution equation: f i x + c i δ t , t + δ t − f i x , t = Ω i x , t + δ t F i , i = 0,1 , … , q − 1 , (11) where f _i ( x , t) is the distribution function for particles at position x and time t with discrete velocity c _i , δ _t is the time step, q is the number of discrete velocities, and Ω _i ( x , t) is the discrete collision operator, and F _i accounts for the body force F .

In LBM, the most widely used collision model is the single-relaxation-time or Bhatnagar-Gross-Krook (BGK) model. However, it has been shown that the BGK model has some shortcomings in pore-scale simulations, such as the unphysical viscosity-dependent permeability. On the other hand, the Multiple-Relaxtion-Time (MRT) model (Heyhat et al., 2020) can effectively solve the problem by introducing different relaxation times (Lallemand and Luo, 2000). Furthermore, the MRT model can also enhance the numerical stability, which is particularly useful for the present study where the viscosity ratio of the two fluids is large. Therefore, we will use the MRT model in the present study.

The collision operator in the MRT model can be expressed as: Ω x , t = − M − 1 S m x , t − m e q x , t , (12) where m and m ^(eq) are the moment and corresponding equilibria in moment space, respectively, and M is a q × q transformation matrix that maps the distribution functions to the moments’ space, m = M ⋅ f , S = diag(s ₀, s ₁, … , s _q−1) is a diagonal matrix of relaxation rates. In this work, we consider two-dimensional problems and use the two-dimensional nine-velocity (D2Q9) model where the discrete velocities are defined by c i = 0,0 , i = 0 , 1,0 c , 0,1 c , − 1,0 c , 0 , − 1 c , i = 1 ∼ 4 , 1,1 c , − 1,1 c , − 1 , − 1 c , 1 , − 1 c , i = 5 ∼ 8 , (13) where c = δ _x /δ _t , with δ _x being lattice spacing. In the present work, c = 1. The transformation matrix M is defined as follows: M = 1 1 1 1 1 1 1 1 1 − 4 − 1 − 1 − 1 − 1 2 2 2 2 4 − 2 − 2 − 2 − 2 1 1 1 1 0 1 0 − 1 0 1 − 1 − 1 1 0 − 2 0 2 0 1 − 1 − 1 1 0 0 1 0 − 1 1 1 − 1 − 1 0 0 − 2 0 2 1 1 − 1 − 1 0 1 − 1 1 − 1 0 0 0 0 0 0 0 0 0 1 − 1 1 − 1 . (14) The corresponding discrete velocity moments of the distribution function are m = ρ , e , ε , j x , q x , j y , q y , p x x , p x y T , (15) where T represents the transpose operator, ρ is the fluid density, e and ɛ are related to the total energy and the energy square, j _x and j _y are components of the momentum, i.e., j _x = ρu _x , j _y = ρu _y , q _x and q _y are the x and y components of the energy flux, p _xx and p _xy are related to the symmetric and traceless components of the stress tensor, respectively. For an incompressible fluid, the density of the fluid is approximately uniform and is denoted by ρ ₀, the density fluctuation is δρ, thus ρ = ρ ₀ + δρ. The corresponding equilibrium expressions of the moments are given by m e q = ρ , − 2 ρ + 3 ρ 0 u 2 , ρ − 3 ρ 0 u 2 , ρ 0 u x , − ρ 0 u x , ρ 0 u y , − ρ 0 u y , ρ 0 u x 2 − ρ 0 u y 2 , ρ 0 u x u y T , (16) and the relaxation matrix corresponding to the nine moments is S = d i a g s ρ , s e , s ε , s j , s q , s j , s q , s ν , s ν . (17) It should be noted that the fluid density and momentum are conserved during the collision process so that the relaxation rates corresponding to these moments, s _ρ and s _j , can take arbitrary values. The other relaxation rates are given by s _e = s _ɛ = s _ν = 1/τ _ν and s _q = 8(2 − s _ν )/(8 − s _ν ) (Pan et al., 2006), where τ _ν is determined by the dimensionless kinematic viscosity ν(C), ν C = 1 3 τ ν − 1 2 δ t . (18)

The body force is defined as follows: F = M − 1 I − 1 2 S M F ̄ , (19) where F = F 0 , F 1 , … , F 8 T , F ̄ = F 0 ̄ , F 1 ̄ , … , F 8 ̄ T , (20) and F i ̄ = ω i c i ⋅ F c s 2 + u F : c i c i − c s 2 I c s 4 . (21)

In the moment space, the body force F ̂ can be derived using Eq. 21 and the transformation matrix M F ̂ = M F = I − 1 2 S M F ̄ , (22) where M F ̄ = 0,6 u ⋅ F , − 6 u ⋅ F , F x , − F x , F y , − F y , 2 u x F x − u y F y , u x F y − u y F x T . (23)

The fluid density and velocity can be obtained through the distribution function: ρ = ∑ i = 0 8 f i , ρ 0 u = ∑ i = 0 8 c i f i + δ t 2 F . (24)

The concentration field is described by the lattice kinetic scheme (Inamuro, 2002), the evolution equation is given as: g i x + c i δ t , t + δ t = g i e q x , t , i = 0,1 , … , 8 , (25) where g _i ( x , t) accounts for the concentration C, and the discrete velocity set is the same as that used in the above MRT model for the velocity field; The equilibrium distribution function g i ( e q ) is defined by g i e q = ω i C 1 + c i ⋅ u c s 2 + ω i A δ t c i ⋅ ∇ C , (26) where the weight coefficients are ω ₀ = 4/9, ω ₁₋₄ = 1/9, ω ₅₋₈ = 1/36, c s 2 = c 2 / 3 , and the parameter A is related to the Schmidt number, 1 S c = c s 2 1 2 − A δ t . (27) The concentration is defined by the distribution function g _i , C x , t = ∑ i = 0 8 g x , t . (28)

The concentration gradient in Eq. 26 can be obtained from the first-order moment of the non-equilibrium function at the given point, − c s 2 δ t ∇ C = ∑ i c i g i − g i e q , (29) and the final result is ∇ C = ∑ i c i g i − C u c s 2 A − 1 δ t . (30)

Implementing boundary conditions is a fundamental problem in LBM (Lou et al., 2018). Since the fluid’s concentration and velocity are known at the inlet of the porous media, the non-equilibrium extrapolation scheme (Guo et al., 2002) is applied. The no-slip boundary condition is realized by the halfway bounce-back scheme (Ladd, 1994). It can be shown that if the relaxation rate τ _q is chosen as s _q = 8(2 − s _ν )/(8 − s _ν ) in the MRT model, the no-slip boundary condition can be realized accurately and spurious slip can be avoided (Pan et al., 2006). For the no-flux boundary condition, a similar bounce-back scheme is used (Wang et al., 2013). For the Neumann boundary condition at the outlet, the corresponding boundary condition scheme for the lattice Boltzmann method is given in Ref Lou et al. (2013).

The above LBM was previously proposed in Ref Liu et al. (2016) for solving Eqs 4–7. It has been validated that the above model has second-order accuracy in space, is insensitive to relaxation parameters, and is very stable at high Péclet number and large viscosity ratio compared with the lattice BGK model. As such, the model can accurately simulate the fluid flow and diffusion in porous media at high Péclet number and large viscosity ratio.

We now use the LBM mentioned above to numerically simulate the miscible displacement in an inclined porous medium at the pore scale. The governing equations and boundary conditions for this problem are described in Sec. 2. The lattice size is 1,600 × 400. The robustness of our results has been tested successfully with longer and wider lattices (from 800 × 200 to 3,200 × 800), while keeping the same characteristic dimensionless numbers. The parameters are set as follows: Ra = 10⁶, Sc = 80, and M = 54.6.

The temporal evolution of a dimensionless measure of the concentration of the displacing fluid “1” is plotted in Figure 2 to start the presentation of our results. The selected parameter values are typical of a situation in which a less viscous fluid displaces a more viscous fluid; in this case, one would expect the flow to be destabilized due to viscous contrasts and porous medium structure. As can be seen in Figure 2, at t = 0.125, the displacing fluid has just entered the porous medium and the fluid-fluid interface is clear. Then, the displacing fluid finds a path of least resistance, and the miscible displacement of fluid 2 by fluid 1 is accompanied by the development of instabilities, these manifest themselves via the formation of fingering structures. At t = 0.750 and 1.375, instabilities of the branching phenomena also arise from the structure of the porous medium. At the latter stages of the flow (t = 2.000), it can be seen that a distinct “dominant band” forms, and fingerings develop along this band in the porous medium.

In this paper, we have studied the displacement process of miscible fluids in porous media at the pore scale while taking gravity into account by using a LBM. The study investigated the temporal distribution of fluid concentration fields in porous media under gravitational influences. Gravity’s impact on the displacement process was explored by comparing it to cases with no gravitational influence. Additionally, the impact of inclination angle (θ) and viscosity ratio (M) on interface stability and displacement efficiency was investigated. The conclusions are as follows.

1) As the inclination angle increases, the viscous fingering instability in porous media becomes more apparent. The displacement efficiency is divided into two regions when the inclination angle increases. One region is when the inclination angle ranges from 0° to 50°. At this time, as the inclination angle increases, so does the displacement efficiency. In the other region, the inclination angle ranges from 50° to 90°, and as the inclination angle increases, the displacement efficiency remains relatively stable.