A lattice Boltzmann exploration of two-phase displacement in 2D porous media under various pressure boundary conditions

Gunxi Yn,Zi Li,Thierry Bore,Sergio Andres Glindo Torres,,Alexnder Scheuermnn,Ling Li,

a School of Civil Engineering,University of Queensland,St Lucia,Brisbane,QLD,4072,Australia

b School of Engineering,Westlake University,Hangzhou,310024,China

c Institute of Advanced Technology,Westlake Institute for Advanced Study,Hangzhou,310024,China

Keywords:Two-phase flow Porous media Dynamic effects Multistep in/outflow Capillary pressure Interfacial area Flow regimes

ABSTRACT While experimental designs developed in recent decades have contributed to research on dynamic nonequilibrium effects in transient two-phase flow in porous media,this problem has been seldom investigated using direct numerical simulation (DNS).Only a few studies have sought to numerically solve Navier-Stokes equations with level-set (LS) or volume-of-fluid (VoF) methods,each of which has constraints in terms of meniscus dynamics for various flow velocities in the control volume(CV)domain.The Shan-Chen multiphase multicomponent lattice Boltzmann method (SC-LBM) has a fundamental mechanism to separate immiscible fluid phases in the density domain without these limitations.Therefore,this study applied it to explore two-phase displacement in a single representative elementary volume (REV) of two-dimensional (2D) porous media.As a continuation of a previous investigation into one-step inflow/outflow in 2D porous media,this work seeks to identify dynamic nonequilibrium effects on capillary pressure-saturation relationship (Pc-S) for quasi-steady-state flow and multistep inflow/outflow under various pressure boundary conditions.The simulation outcomes show that Pc,S and specific interfacial area (anw) had multistep-wise dynamic effects corresponding to the multistep-wise pressure boundary conditions.With finer adjustments to the increase in pressure over more steps,dynamic nonequilibrium effects were significantly alleviated and even finally disappeared to achieve quasisteady-state inflow/outflow conditions.Furthermore,triangular wave-formed pressure boundary conditions were applied in different periods to investigate dynamic nonequilibrium effects for hysteretical Pc-S.The results showed overshoot and undershoot of Pc to S in loops of the nonequilibrium hysteresis.In addition,the flow regimes of multistep-wise dynamic effects were analyzed in terms of Reynolds and capillary numbers(Re and Ca).The analysis of REV-scale flow regimes showed higher Re(1 ＜Re ＜10)for more significant dynamic nonequilibrium effects.This indicates that inertia is critical for transient twophase flow in porous media under dynamic nonequilibrium conditions.

1.Introduction

The theory of multiphase flow in porous media has applications for the transport of nonaqueous-phase contamination in environmental hydrogeology,the failure of unsaturated soil slopes in geotechnical engineering,and the enhancement of oil recovery(EOR) in petroleum engineering.Analytically and numerically solving this theory can help predict the migration of contamination in unsaturated soil and the variations in shear strength of unsaturated soil with soil moisture dynamics as well as the oil that can be recovered by injecting another immiscible fluid into a deep reservoir.It is,therefore,a crucial geophysical process in need of comprehensive investigation.The conventional theory is an extension of the groundwater equation by counting the relative permeability and specific storage in the function of either the capillary pressure (Pc) or saturation (S) (Bear,1972).Mathematically,the partial differential equation in diffusion form is extended to nonlinear diffusion form with a diffusivity term in the function of the state variables(Richards,1931;Bear,1972).Thus,the prediction of both S and Pcdepends on iteratively calculating the state variables on each representative elementary volume(REV)by applying experimentally determined constitutive relationships (Pc-S,relative permeability-saturation Kr-S) and suitable hydraulic boundary conditions specified for replicating the system in practice.The constitutive relationships are the foundations of the theory because both state variables of this pair are determined by Pc-S.Moreover,with the development of pedo-transfer functions between Pc-S and Kr-S (Brooks and Corey,1964; Mualem,1976,1978; Van Genuchten,1980; Fredlund et al.,1994),Kr-S,as a curve demanding complexity in experimental testing,can be quickly determined from Pc-S curve.Hence,determining Pc-S relationship is the critical core for preserving the accuracy of modeling of twophase flow seepage.

The Pc-S relationship is typically measured using the axis translation technique(ATT),which involves the displacement of two immiscible-phase fluids in a pressure cell with a high-air-entry(HAE) ceramic plate under a specimen to prevent gas leakage(Fredlund and Rahardjo,1993;ASTM D6836-02,2002;Lu and Likos,2004).According to the testing standards for unsaturated soil(ASTM D6836-02,2002;ASTM D7664-10,2010),the hydraulic resistivity of this plate is usually several magnitudes greater than that of the soil specimen.Setting up a boundary with such low permeability makes it less likely that the Pc-S relation can be recorded under the conditions of nonequilibrium transient flow.The physical process for such flow conditions can be described as a high-pressure fluid in the nonwetting phase pushing a positively pressurized wetting-phase fluid out of the flow cell without hydraulic energy reduction on the boundaries.The dynamic Pcor nonequilibrium soil suction,therefore,has to be monitored by tensiometers.The dynamic nonequilibrium Pc-S relationwas discovered 50 years ago in a series of pioneering investigations (Topp et al.,1967; Stauffer,1978; De Gennes,1988).In these studies,high-velocity two-phase seepage was generated by increasing boundary flow rates.As a result,these researchers determined that the dynamic nonequilibrium Pc-S relationship was dependent on velocity (Kalaydjian,1992;Wildenschild et al.,2001;Hassanizadeh et al.,2002;O’Carroll et al.,2005;Schembre and Kovscek,2006;Mirzaei and Das,2007;Sakaki et al.,2010; Das and Mirzaei,2012; Yan et al.,2021a; Yan et al.,2022a).Specifically,they found that deviations between the dynamic nonequilibrium and static/steady-state Pc-S measurements could be enlarged at high seepage velocities.Experiments on twophase displacement were also used to validate the Richards’equation(Richards,1931)and the theory of two-phase displacement in porous media (Bear,1972),which failed to simulate these experimental outcomes (Liakopoulos,1964; Hassanizadeh et al.,2002).The physical reasons for this were discussed but had remained uncertain by that moment.This subsequently led to modern research on dynamic multiphase flow seepage from the pore scale up to single-REV/continuum scales(Joekar-Niasar et al.,2010a,b;Joekar-Niasar and Hassanizadeh,2012a,b; Karadimitriou and Hassanizadeh,2012;Karadimitriou et al.,2014;Ferrari et al.,2015).

Since its discovery,several theories have been developed regarding the dynamic Pc-S.Barenblatt and Vinnichenko (1980)developed dynamic moisture theory to describe the nonequilibrium in the form of moisture/saturation relaxation.The mathematical formulation empirically elucidated the moisture/saturation dynamics under transient flow conditions but neglected variations in the nonequilibrium capillary pressure from a physical perspective.Hassanizadeh and Gray (1993b) and Kalaydjian (1987) subsequently provided a thermodynamic-based theory that considers the free energy released from the interfaces of immiscible fluids.However,this complex theory includes material coefficients timed by the saturation and specific interfacial area as additional driven terms in a Darcian form of momentum balance,and an extra production term as a source-sink term in the continuity equation of specific interfaces(anw).These terms led to significant challenges with regard to the experimental instruments used (Hassanizadeh and Gray,1993a,b).To simplify the theory,Hassanizadeh et al.(2002) developed an equation of dynamic capillary pressure at REV scale to characterize the relaxation in the nonequilibrium capillary pressure with variations in saturation by adding only a single dynamic coefficient that was experimentally detectable.The inclusion of this equation resulted in a third-order term when applied to the conventional theory (Hassanizadeh et al.,2002).Since then,most soil column and pressure cell tests have validated the dynamic coefficient as well as its dependence on the flow velocity,size of the specimen,temperature,permeability,pore geometry,etc.(O’Carroll et al.,2005; Mirzaei and Das,2007; Sakaki et al.,2010; Das and Mirzaei,2012; Hanspal and Das,2012;Abidoye and Das,2014; Zhuang,2015; Zhuang et al.,2016,2017,2019).Comprehensive reviews of prevalent issues in dynamic nonequilibrium in the context of theories,experiments,and macro/microscale mechanisms have been provided by Diamantopoulos and Durner (2012),Durner et al.(2014),and Yan et al.(2022b).

The above studies and reviews have focused merely on continuum-scale experiments and modeling.Several reasons for the dynamic nonequilibrium effects have been proposed and discussed without investigating them at micro/mesoscale due to observational constraints in conventional testing and modeling techniques.To overcome these constraints,some methods have been developed at pore scale and REV scale in the last two decades,including pore network modeling(Joekar-Niasar et al.,2008,2009,2010; Joekar-Niasar and Hassanizadeh,2012a,b; Sweijen et al.,2018; Hu et al.,2019; Huang et al.,2020),physical microscale models (Karadimitriou and Hassanizadeh,2012; Karadimitriou et al.,2012,2013,2014; Kunz et al.,2015),volume-of-fluid (VoF)-based direct numerical simulations(DNS)(Ferrari and Lunati,2013,2014; Ferrari,2014; Ferrari et al.,2015; Palakurthi et al.,2018;Konangi et al.,2021),level-set (LS)-based DNS methods (Helland et al.,2017; Wang et al.,2017),and multiphase smooth particle hydrodynamics (SPH)-based DNS methods (Sivanesapillai and Steeb,2016,2018; Sivanesapillai et al.,2016).In addition to them,the multiphase lattice Boltzmann method (LBM) is a powerful numerical tool.Shan-Chen multiphase multicomponent LBM(SCLBM)applies three pseudo-interactions to the fluid density domain to reproduce the fluid phase separations (Shan and Chen,1993;Sukop and Thorne,2007).Thus,SC-LBM has intrinsic advantages for investigating dynamic nonequilibrium two-phase flow in porous media.SC-LBM has been applied to study Pc-S,Kr-S and Pc-S-anwconstitutive relationships under static and steady-state conditions in a series of numerical studies (Pan et al.,2004;Schaap et al.,2007;Ahrenholz et al.,2008;Liu et al.,2015;Galindo-Torres et al.,2016;Li et al.,2018,2019).This method has also been validated by physical experiments by Porter et al.(2009).

However,SC-LBM has been rarely applied to investigating dynamic nonequilibrium effects in comparison with other pore-scale investigations.The earliest work on ideal two-dimensional (2D)porous media was from Porter et al.(2006),who simulated the dynamic nonequilibrium Pcand Pc-S under transient flow conditions using SC-LBM.The outcomes showed many fluctuations of Pcin the dynamic Pc-S only for primary drainage.Galindo-Torres et al.(2013) applied the single-component multiphase SC-LBM to investigate dynamic effects in the hysteretical water retention behavior in a 2D porous media.The simulation outputs confirmed the so-defined"hydraulic ratcheting"phenomenon experimentally discovered by Scheuermann et al.(2005),and later verified by Scheuermann et al.(2014).However,with an aim to study groundwater dynamics by counting the saturated capillary fringe inside the groundwater table,the water potentials were represented in the elevation head instead of by Pcor soil suction.Yan et al.(2017) continued this work using a multicomponent multiphase SC-LBM to simulate dynamic effects in Pc-S-anw.The results confirmed the previous experimental finding whereby anwcould be enlarged in highly nonequilibrium conditions for both drainage and imbibition.Later,Yan et al.(2018) further conducted an SC-LBM simulation by applying various boundary conditions,and identified the dependence of dynamic Pc-S on pressure boundary conditions.Meanwhile,Tang et al.(2018) investigated dynamic effects in Pcfor two-phase fracture flow,followed by an examination of the dynamic Pcin sandstone using SC-LBM (Tang et al.,2019; Cao et al.,2020).Although their study focused on natural geomaterials,the dynamic Pcwere reproduced using SCLBM,followed by a determination of dynamic coefficients.Nevertheless,owing to the high computational cost of SC-LBM,their studies still concentrated on primary drainage.Recently,Yan et al.(2021b) investigated dynamic nonequilibrium two-phase inflow/outflow in quasi-2D porous media for both drainage and imbibition by applying one-step-imposed pressure boundary conditions.The most comprehensive reviews of pore-scale modeling of dynamic effects have been provided by Chen et al.(2022) and Yan et al.(2022b).However,the nonequilibrium two-phase multistep inflow/outflow has not been comprehensively explored using SCLBM by far.

Therefore,this study investigates nonequilibrium transient twophase flow in porous media using SC-LBM.The multistep and triangular wave-formed pressure boundary conditions were applied to generating dynamic Pc-S in the forms of multistep-wise and hysteretical loops.In addition,the dynamic multistep-wise Pc-S against anwand seepage flux (q) as calculated by the simulation are discussed by referring to the neglected inertial terms in the multiphase Navier-Stokes (NS) equations.Compared with previous experimental studies on dynamic nonequilibrium effects from a velocity-dependent perspective based on thermodynamic-based theories,this study provides a distinctive vision for exploring this issue from a pressure-dependent perspective with a multiphase hydrodynamic basis.In principle,this numerical investigation advances the understanding of dynamic multiphase inflow/outflow through porous media at a single-REV scale and micro/mesoscale.

2.Numerical methodology

2.1.Shan-Chen LBM

The LBM is a computational method that numerically approximates the solution of the Boltzmann equation within a discretized domain of the particle probability function and corresponding velocity.There are several ways to discretize the simulation domain.Here,the D2Q9 scheme was selected because this study investigates only dynamic two-phase displacement in porous media as a 2D problem.The discretized lattice Boltzmann equation shows that a streaming process is the consequence of a collision process whereby particles arrive in and out of a target region (Sukop and Thorne,2007),as follows:

where fiis a function of the probability density for finding molecules at location x at node i,t is the time,eiis the velocity of the discrete fluid at node i,Δt is the time step,fieqis the Maxwell-Boltzmann equilibrium distribution function,and τ is the characteristic relaxation time.The D2Q9 discretizing scheme for eiis shown in Fig.1 (Sukop and Thorne,2007).

Fig.1.The D2Q9 cell depicting the discrete fluid velocities (ei,i=0-8) for a located probability function fi.

To solve Eq.(1) in a domain consisting of the cells shown in Fig.1,the macroscopic fluid density ρ and macroscopic fluid velocity v for each node i in the streaming process on the left-hand side are given by

and the discrete velocity(ei)is calculated by the following matrix:

For the collision process on the right-hand side in Eq.(1),fieqis derived from the Taylor series of the Maxwell-Boltzmann distribution function to recover the NS equations of the hydrodynamic model from the Boltzmann equations (Sukop and Thorne,2007).Thus,the Bhatnagar-Gross-Krook(BGK)approximation is applied as

where veqis the equilibrium velocity;Csis the speed of sound in the lattice,and Cs=(the basic speed in lattice is C=Δx/Δt,the lattice unit length Δx=1 lattice unit (lu) and the lattice unit time step Δt=1 time step(ts));and the weighing factor(wi)can be given by

The relaxation time τ determines the numerical stability of LBM.The kinematic viscosity (υ) is a function of τ:

where in order to maintain numerical stability,τ should be much higher than 0.5 to deliver a positive υ with a physical meaning.Based on work by Sukop and Thorne (2007),numerical instability also occurs when τ is close to 0.5.Hence it is recommended to be close to one.The υ was thus set to 0.16 lu2/ts.The single relaxation time(SRT)was selected instead of the multi-relaxation time(MRT)because this work did not consider the effects of high density and viscosity ratios on dynamic nonequilibrium.Applying an SRT with unit density and viscosity ratios ensures numerical stability and model validity.In addition,the veqis calculated by

where F is the total external force.For the single-phase singlecomponent LBM,F is the gravitational force ρg for simulating problems on the continuum scale(more than one REV).This study simulated two-phase flow for a single REV of porous media,because the gravitational acceleration g was set to zero.Using Eqs.(1)-(8),the hydrodynamic model can be numerically solved by the single-phase LBM.

In addition to single-phase hydrodynamics,Shan and Chen(1993)introduced three pseudo-interactions and multiple components.Therefore,the SC-LBM can be used to separate a fluid from another by setting interactions between each pair of fluid phases (fluid-fluid,fluid-solid).Each fluid phase has to be allocated with independent density and velocity fields.Moreover,nucleation and dissolution can be replicated by this method.Fluid-fluid phase separation is mathematically yielded by the pseudo-repulsive interaction:

where Fris the repulsive interaction between the wetting and nonwetting phases,Gris the repulsive intensity controlling the two-phase repulsive interaction,ρj,?is the density of a fluid phase(?),and it is expelled out by the other fluid phase(?′)with density ρj,?′.To generate the hydrophilicity and hydrophobicity for the wetting and nonwetting phases,respectively,the interaction between each fluid and solid phase is calculated by

where Fsis the interaction between a fluid (? or ?′) and a solid; a positive Gsdenotes the hydrophobicity of the nonwetting phase,a negative Gsdenotes the hydrophilicity of the wetting phase; s is a logic function that highlights adjacent solid nodes,and it is switched to zero if no solid nodes are in contact with fluid nodes.Having included these additional external interactions in Eq.(8),F should be redetermined as

Furthermore,due to the external forces and multiple components in Eq.(8),the macroscopic velocity v of the two fluid phases should be corrected by

Thus,with the addition of Eqs.(9)-(12)and two fluid phases in the simulation domain,immiscible multiphase separation and mixing behavior can be simulated by SC-LBM.

2.2.Upscaling from lattice mesoscale to macroscale

The LBM simulations were carried out using the open-source library Mechsys to solve the multiphase multicomponent SC-LBM(http://mechsys.nongnu.org,accessed in 2021).This numerical package is still under development by Prof.Sergio Andres Galindo Torres at the School of Engineering at Westlake University in Hangzhou,China,and at the Geotechnical Engineering Centre(GEC) of the University of Queensland (UQ) in St.Lucia,QLD,Australia.This study applied the D2Q9 Shan-Chen multiphase multicomponent SC-LBM simulation.The open-source LBM solver can output the density and velocity fields for each time step.Many previous studies have verified SC-LBM on the Mechsys open-source library (Galindo-Torres et al.,2013,2016; Li et al.,2018,2019; Yan et al.,2021b).Therefore,this library can be considered a mature SC-LBM model for geoscientific applications.As the simulations were implemented on a single REV of 2D granular porous media,the macroscale state variables needed to be determined by the following methods:

(1) The pressure(Pj,?)of a fluid phase ? on node j was calculated by the equation of state (EOS) of SC-LBM:

(2) Then,the pressures of both fluid phases on the entire simulation domain could be upscaled by

where Aj,? is the volume at node j occupied by phase ?.It denoted an area for the setup of the 2D simulation.According to the wettability settings,the wetting and nonwetting phase could be assigned to either phase ? or ?′.

(3) The capillary pressure(Pc) was then given by

where n and w indicate the nonwetting and wetting phases,respectively.

(4) The saturation (S?) of each fluid phase is

and,based on such a density threshold,any density value in any other range can be identified as a combination of interfacial area and the other fluid phase while excluding the solid phase(ρj=0).

(5) The specific interfacial area is

which is the interfacial area (Anw) over the entire 2D simulation domain (AT).The isosurface (·) function in MATLAB (MathWorks Inc.,Natick,MA,USA)was applied to detecting and calculating the two-phase interface in the quasi-2D domain.

(6) The velocity field was upscaled to Darcy flux (q) using

where n is the porosity and vjis the velocity on node j.

(7) For dynamic multiphase flow,the bound number (Bo),capillary number(Ca),and Reynolds number(Re)are usually considered to analyze the flow regimes for a single REV of porous media.As this simulation focused on a single REV,gravity was not included for the zero-dimensional scenario,and Bo did not need to be considered.Re and Ca were separately calculated by

where deff,satis the effective hydraulic character length for singlephase seepage under fully saturated conditions,and σ is the surface tension of the wetting phase.The effective hydraulic character length for multiphase seepage needed to decrease with the history of saturation S(t).The geometric mean of the particle size distribution was used in this 2D simulation.

(8) The velocities on the two-phase interface near the solid boundaries and particle surfaces are unreliable because of spurious flow velocities (Chen et al.,2014,2018).Thus,a filtering process was needed (Pooley and Furtado,2008).A marching algorithm was then applied in MATLAB to eliminating the densities and velocities in the first layer of adjacent nodes around the solid particles.

2.3.Simulation setup

Numerically approximating the solution of SC-LBM is a computationally demanding job.Even though Mechsys was used to solve it in the 2D domain,a single simulation usually took several days on the authors’ desktop.It could be completed in two days using the high-performance supercluster Goliath,available at the Faculty of Engineering of Architecture and Information Technology(EAIT)of UQ.The fluid density and viscosity ratios were set one to reduce the computational cost by calculating through

where Pcphyiscis the capillary pressure of a physical fluid,Pc,Lis the lattice capillary pressure,σphysicis the surface tension of the physical fluid,σLis the lattice surface tension,ρphysicis the density of the physical fluid,μphysicis its dynamic viscosity,υ is the lattice kinematic viscosity,Δxphysicis the physical unit length,Δtphysicis the physical unit time step,and Δt is the lattice unit time step.The simulation outcomes can be then presented according to the lattice-physical unit transfers in Table 1.

Table 1 The lattice-physical unit transfers.

As the selection of Grand ρjdetermined the repulsive interaction between two immiscible fluids and the compressibility of each fluid,the contact angle was set to zero,to simplify the problem and alleviate the computational cost,by

With this setting,the lattice surface tension can be given by the Young-Laplace equation:

where Pinand Poutare the pressures of fluid components inside and outside a single simulated bubble,respectively; and the surface tension σ is given in Tables 1 and 2.Same fluid properties were used by Galindo-Torres et al.(2016),Yan et al.(2017,2021b),and Li et al.(2018,2019).In addition,the details of choosing the density and pseudo-force have been provided by Shan and Chen(1993),Martys and Douglas (2001),and Schaap et al.(2007).

The quasi-2D simulation domain had a size of 5000 μm 5000 μm 10 μm(500 lu 500 lu 1 lu),and was filled with a package of granular disks.Based on the criteria for the mesh from Sukop and Thorne (2007) and Pan et al.(2004),at least 5 lu between two solid boundaries were required to ensure the validity of flow in a simulation,and the simulation domain needed to be independent of its mesh size (lx=ly) when l was 13 times greater than the mean grain size(D50).This porous package met both criteria because of the pore channels over 5 lu and l/D50=5000 μm/360 μm=13.89.Fig.2 shows that to avoid the loose compaction of the porous media that leads to local heterogeneity,a few particles near the boundaries were cut by the edges of the domain.The image analysis algorithm developed by Rabbani et al.(2014)was applied to measuring the grain size distribution (GSD).The single-phase LBM was also carried out for the same 2D particle package to determine the intrinsic permeability Ksat.A summary of specifications of the porous media and fluid properties is listed in Table 2.This simulation setup can be considered to model two-phase flow displacement through an unconsolidated REV of uniform medium sand,although the natural sand particles are not perfectly spherical.

Table 2 Parameters settings for the numerical simulation.

This simulation aimed to replicate two-phase transient displacement as a zero-dimensional point-wise REV in a column of porous media.In comparison with conventional pressure cells,it had no ceramic material with low permeability under the specimen.Two solid boundaries on the left and right were simulated using the bounce-back scheme (Sukop and Thorne,2007).Two functions were applied to assigning the densities at the top and bottom of the domain to simulate quasi-steady-state and transient flow conditions.With a density ratio of one,the density added at the top for the nonwetting phase was equal to the density reduced at the bottom for the wetting phase in order to conserve the fluid mass for simulating drainage,and vice versa for simulating imbibition.To simulate multistep inflow/outflow,the density function is given by

where ρ?,iis the density of phase component (?) assigned to time step i;ρ?,0is the initial density;Δρ is the incremental step in density for each time step;each time step is given by taking the floor of the ratio of the current simulation time(ti)to the incremental time step Δt,which is the equal to the total time(Tf)divided by the number of increments in density (Nstep).

To model the quasi-steady-state inflow/outflow,the infinitesimal increment in density is given by the triangular wave function:

where a=π/ω and ω=2πNcycle/Tf,in which Ncycleis the number of hydraulic cycles.Because two components were used to simulate phase separation and mixing,the density control for each fluid phase at the top and bottom boundaries should follow

where ρf,bis the density assigned to each fluid phase f(f=w/nw for wetting or nonwetting)on the boundary b(b=top/bottom for the top or bottom boundary).As the increasing and decreasing densities can be transferred to the increasing and decreasing pressures using LBM EOS (Eq.(13)),these two density boundary conditions can be seen as Dirichlet boundary conditions for a hydrodynamics simulation.Assigning densities to the top and bottom boundaries to replicate pressure boundary conditions has been applied by Pan et al.(2004),Schaap et al.(2007),Porter et al.(2009) and Galindo-Torres et al.(2016).By following those previous experiences,the transient and steady-state flow conditions could be reproduced to investigate the responses of the state variables and the constitutive relationships in the theory of multiphase seepage flow.The multistep and triangular wave functions of density are respectively shown in Fig.3a and b.

3.Results and discussion

3.1.Dynamic Pc-S curves under transient multistep inflow/outflow

To examine dynamic multiphase flow seepage in porous media,both the capillary pressure(Pc)and saturation of the wetting phase(S) against time (t) are plotted in Fig.4.The transient flow conditions were generated by increasing or reducing pressure on the boundaries in several steps.As shown in Fig.4a and d,the curves of S-t for both drainage and imbibition in a few steps of pressurization (1-5 steps) were significantly distinct from those obtained in many more steps (10-50 steps).Moreover,with the increasing number of steps of pressurization and decreasing increments in pressure in each step,the curves of S-t finally overlapped,indicating that such a slow multistep inflow/outflow tended to satisfy the assumption of instantaneous equilibrium.Furthermore,for the Pc-t in Fig.4b and e,the staircase-shaped curves similar to the step-wise density functions in Fig.3a were observed for multistep drainage and imbibition.The abruption of Pcwas also noted for one-and two-step inflows/outflows when higher pressure increments were applied on the boundaries.Such a Pcvariation corresponded to the variation of pressure boundary conditions for multistep inflow/outflow illustrated in Fig.3a.Regarding this illustration,the dynamic Pc-t might have been partially dependent on the variation of pressure boundary conditions.The nonequilibrium transient flow could be simulated with fewer steps to vary the pressure and higher increments in it.Such dynamic effects could be alleviated by increasing the number of steps and reducing the incremental interval of pressure in each step.In summary,the smaller the increments in pressure were in each step when more steps were used (10-50),the greater the extent to which the Pc/S-t datasets matched each other was,and the more the quasi-steadystate flow conditions approached.

The dynamic Pc-S curves are plotted in Fig.4c and f for both drainage and imbibition based on the curves of Pc-t and S-t.Like many previous experimental and numerical findings (Topp et al.,1967; Hassanizadeh et al.,2002; O’Carroll et al.,2005; Schembre and Kovscek,2006; Das and Mirzaei,2012; Yan et al.,2021b),where the dynamic Pcovershoot for drainage and undershoot for imbibition,these multistep inflow/outflow simulations confirmed this difference between static/quasi-steady-state and dynamic Pc-S curves.For one-and two-step pressure increments corresponding to each value of S,the dynamic Pcwas larger than Pcfor steady-state drainage and smaller than Pcfor steady-state imbibition.Because the pressure boundary conditions were varied in a step-wise manner,the deviations between the steady-state and dynamic Pc-S showed similar step-wise behavior.When the number of steps of pressurization was increased from 10 to 50,the dynamic Pc-S curves were closer to those generated under the steady-state flow conditions.These results also confirm the experimental findings by Topp et al.(1967).The dynamic curves reverted to steady-state/static curves when the seepage flow velocity decreased to achieve steady-state flow conditions.As shown in Fig.4c and f,the equilibrium condition was reached after the pressure on the boundaries was increased for a certain period.Once the pressure rose again,the second transient flow was reproduced such that a deviation between the equilibrium and nonequilibrium responses of Pcwas generated again.However,the static/steady-state curves determined by the conventional testing method entirely rely on selecting the data points for the static Pcand S or generating quasisteady-state flow with a very low velocity.Whether such a steadystate flow and instantaneous equilibrium conditions can be conserved highly relies on the dynamic boundary conditions as a function of time.According to this multistep inflow/outflow simulation,the dynamic effects in Pc-S curves could be a function of dynamic pressure boundary conditions.Moreover,their uniqueness under transient multistep inflow/outflow can be ensured only if the pressure boundary conditions are adjusted slowly with an infinitesimal variation in pressure in each step.

Fig.2.SC-LBM simulation setup:(a)The quasi-2D granular porous media used in the numerical simulations,and(b)The grain size distribution of the quasi-2D porous media in the simulation domain.

Fig.3.The demonstration of the pressure boundary conditions applied by assigning the variations in density (ρ?,i) at the top and bottom boundaries: (a) Multistep incremental function of density(Eq.(25)),and(b)Triangular wave function of density(Eq.(26)).A longer total time(Tf)yields a smaller increase in density to ensure instantaneous equilibrium under quasi-steady-state inflow/outflow,while a shorter Tf yields dynamic effects under transient inflow/outflow.

3.2.Flow pattern differences between steady-state and transient inflow/outflow

To explore the physical reasons for differences between the dynamic and steady-state Pc-S curves,a few snapshots of wetting phase fluids in porous media given by the two-step and fifty-step outflow simulations are shown in Table 3 for S=86%,43% and 28%.Two different fluid distributions between the quasi-steadystate (50 steps) and multistep transient (2 steps) flow were obtained for each S.When the wetting fluid was drained out to S=86%,there were fewer distinctions between the two cases.The differences were evident when the fluid was drained down to S=43% and 28%.Compared with the drainage in 50 steps,which showed more isolated trapping of the wetting phase,the drainage in only two steps left more connected fluid phases at the center of the quasi-2D porous media.However,the preferential flow paths of the nonwetting phase exhibited similar flow patterns between the steady-state and transient flow conditions.This series of snapshots shows that the flow patterns for the quasi-steady-state and transient outflows were not precisely the same.For steady-state drainage,the wetting phase fluid preferred to be isolated and trapped inside the pore matrix.By contrast,the wetting phase fluid was usually connected in a region for transient drainage.Such different fluid distributions might have been a physical reason for the overshooting Pcin dynamic drainage curves because of a different value of the localized Pc.These microscale/mesoscale findings agree with many previous multistep outflow tests(Poulovassilis,1974; Wildenschild et al.,2001; Weller et al.,2011).The relevant studies have claimed that for transient drainage,a higher dynamic Pcfor each value of S occurred due to the entrapment of water and its blockage through the pores(Diamantopoulos and Durner,2012).In fact,such a difference in fluid distributions between the transient and steady-state flow conditions was related to the fluid properties and particle geometry.The fluid properties in Table 2 show that the contact angles were fixed to zero without changing wettability during the simulation of outflow.The dynamic interface might have made a minor contribution to the contact angle and interface reconfigurations.Therefore,such entrapment of water and distinct flow paths should depend more on local heterogeneity (particle/pore geometry) at a length scale below a single REV scale (Diamantopoulos and Durner,2012).Similar findings were recently presented in a few numerical investigations of dynamic effects induced by pore geometry or microscale heterogeneities (Mirzaei and Das,2007; Joekar-Niasar et al.,2010a,b).

Fig.4.Dynamic capillary pressure (Pc) and saturation (S) for multistep inflow/outflow simulation: (a) Dynamic S-t for drainage,(b) Dynamic Pc-t for drainage,(c) Dynamic Pc-S curves for drainage,(d) Dynamic S-t for imbibition,(e) Dynamic Pc-t for imbibition,and (f) Dynamic Pc-S curves for imbibition.

Moreover,as the pressure boundary conditions were applied by varying the density along the boundaries,the different density distributions of the fluid inside the pore matrix might have been the other physical reason for a higher Pcfor each value of S.For instance,portions of the wetting phase fluid near the interfaces are highlighted in lighter yellow for two-step transient drainage at S=86% and 28% in Table 3.The density variations around the interface indicated the two-phase mixing,separation and compression of each fluid phase adjacent to the interface.Because SC-LBM is based on solving the momentum balance at a mesoscale between the control volume (CV) and a single fluid molecule,it cannot ensure that there is no localized difference in the fluid densities at each node.Therefore,such a spatial variation in density near the interfaces might have contributed to the dynamic Pc.However,whether this is a physical fact or a numerical artifact cannot yet be experimentally validated.Gray and Hassanizadeh(1991) have theoretically challenged the density of soil water as a constant.Thus,for multistep transient drainage,the overshooting Pcfor each value of S might also have occurred due to differences in the localized fluid density.

Table 3 Comparison of fluid distributions of wetting phase between fifty-and two-step outflow simulations.The wetting phase is shaded in yellow,and the nonwetting phase is shaded in cyan.

The snapshots of wetting phase fluids in porous media in twostep and fifty-step inflow simulations for S=19%,58% and 83%are shown in Table 4.The differences between the two simulations were insignificant for S=19%.When S approached 58%,it is evident that the quasi-steady-state(50 steps)inflow generated a relatively sharp wetting front,whereas finger flow occurred for transient (2 steps) inflow.Once the wetting phase fluid had been injected to achieve S=83%,the transient inflow almost reached final equilibrium.Thus,there were minor distinctions between the cases of quasi-steady-state and transient two-step inflow imbibition.

Table 4 Comparison of fluid distributions of wetting phase between fifty-and two-step inflow simulations.The wetting phase is shaded in yellow,and the nonwetting phase is shaded in cyan.

The convex advancing interfaces with negative contact angles were observed for the transient imbibition.By contrast,the concave advancing interfaces with positive contact angles were noted for the quasi-steady-state imbibition.The concave interfaces usually appear in the pore throats of a porous media because the conventional understanding of two-phase flow is based on the Young-Laplace equation for capillarity.In addition,the dependence of dynamic contact angles on the flow velocity has been theoretically studied in the capillary tube (Sheng and Zhou,1992).However,to the best of the authors’knowledge,the convex interfaces have been rarely observed in previous experimental works.Such an interface implies that the flow predominantly occurred not only due to capillary forces or surface tension,but also due to inertia(Liakopoulos,1964; Raats,1965,1972; Raats and Klute,1968;Sposito,1980;Yan et al.,2021b;Chen et al.,2022).When the mass of the wetting phase fluid at high pressure was injected into the single REV in a step-wise manner,the advancing interfaces received inertia transferred through the entire fluid domain.Therefore,the concave interfaces were reversed to convex.

Not only were the interface dynamics identified,the flow patterns of quasi-steady-state and transient flows also exhibited significant differences for S=58% in Table 4.The capillary-dominant steady-state flow ensured quasi-static interfaces and preferentially seeped through the pore throats with small diameters.However,when inertia became predominant,the wetting phase fluid invaded pores with larger throats more quickly than those with smaller ones.Hence,the finger-flow pattern was observed for transient flow conditions while the sharp wetting front was better conserved for quasi-steady-state flow conditions.In summary,the different flow patterns occurred between these two flow conditions due to a combination of inertia (Chen et al.,2022) and local heterogeneity (particle/pore geometry) (Diamantopoulos and Durner,2012).

Moreover,as for transient flow (2 steps) in Table 4,the fading and darker colors on the fingertips of the wetting phase fluid exhibited localized density variations around the interfaces.As aforementioned,SC-LBM simulates two-phase flow within the density and velocity domains.The density variations at local nodes are inevitable,particularly around the two-phase interfaces,as either a physical fact or numerical artifact.The constancy of water density,in this case,cannot be experimentally verified at present.This might be part of the reason for the dynamic Pc,but still debatable.

Based on above analysis for Tables 3 and 4 as well as Fig.4,it is straightforward to find that the dynamic nonequilibrium effects under transient flow conditions,which were manifested in forms of different Pc,distinctive fluid patterns,interface dynamics and local density variation near interfaces,occurred between a prior and posterior equilibrium states.The differences decreased when the given state was close to the previous or an upcoming equilibrium state.This microscale/mesoscale observation shows that the dynamic nonequilibrium effects might have been related to local heterogeneity (particle/pore geometry),interface dynamics,and variations in the density around the interfaces as well as variations of pressure boundary conditions.

Fig.5.Dynamic capillary pressure (Pc) and saturation (S) for cyclic transient and steady-state inflow/outflow when the total periods of the simulation Tf=4000 ms,800 ms and 8 ms:(a)Dynamic and steady-state Pc-t,(b)Dynamic and steady-state St,and (c) Dynamic Pc-S curves under cyclic transient flow (dynamic hysteresis,Tf=800 ms and 8 ms) in comparison with the steady-state flow (Tf=4000 ms).

3.3.Comparison between steady-state and transient cyclic inflows/outflows

In addition to the multistep inflow/outflow simulations,cyclic inflow/outflow simulations were implemented using the triangular-wave density function.As the instances of boundary conditions shown in Fig.3b,the total periods of the simulations(Tf)were set to 4000 ms,800 ms and 8 ms.A comparison of the S-t,Pc-t and Pc-S curves for these three periods is shown in Fig.5.Fig.5b shows that as Tfwas reduced from 4000 ms to 8 ms,the corresponding period for a single triangular wave decreased.Such a cyclic loop of S-t could not cover 0%-100%of S.As shown in Fig.5b,the smaller the cyclic period was,the smaller the range of S that could be simulated.On the contrary,the Pc-t curves for all three periods in Fig.5a show that the complete range of Pccould be covered.As a result,when both Pcand S were correlated and plotted in Fig.5c,the dynamic hysteretical Pc-S curves under transient flow conditions (Tf=800 ms and 8 ms) also showed a higher Pcfor hysteretical drainage and a lower Pcfor hysteretical imbibition than the possible static/steady-state hysteretical scanning loops enclosed within the steady-state Pc-S curves (Tf=4000 ms).

The above results are the first simulations of dynamic hysteretical scanning loops using a CFD package for two-phase flow.The dynamic hysteretical Pc-S curves for the transient cyclic inflow/outflow test were generated by a cyclic triangular-wave density function with different values of Tfat the same amplitudes(Δρ),as illustrated in Fig.3b.The boundary conditions of dynamic pressure led to transient flow with a higher seepage velocity or variations in S.For multistep and cyclic inflow/outflow tests,the dynamic Pc-S curves deviated from the static/steady-state curves for drainage,imbibition and hysteretical scanning.According to the literature(Wildenschild et al.,2001;Hassanizadeh et al.,2002),the dynamic curves are dependent on the seepage velocity or variations in S.Nevertheless,as they were generated by dynamic boundary conditions,those curves should have been originally dependent on the dynamic pressurizing functions applied to the boundary conditions instead of the seepage velocity inside the testing system,thereby simplifying the causes of the dynamic effects.Given that the uniqueness of the constitutive relationship (Pc-S) depends on the boundary conditions,a simulation that uses the conventional theory of dynamic multiphase flow seepage relies entirely on quasisteady-state boundary conditions to reassemble a series of frames of steady-state flow.If the boundary conditions are varied more quickly to violate the assumption of instantaneous equilibrium(Liakopoulos,1964; Yan et al.,2021b),the Pc-S curves become nonunique.Thus,the Pc-t and S-t curves cannot be simulated by the conventional theory of multiphase flow in porous media.

3.4.Dynamic Pc-S-anw curves under transient multistep inflow/outflow

Beyond the dynamic Pc-S curves,the specific interfacial area(anw)was also considered because Hassanizadeh and Gray(1993b)introduced it as a state variable to show a unique Pc-S-anwconstitutive surface enclosing all possible hysteretical loops.Such uniqueness of the constitutive surfaces was confirmed by Porter et al.(2009) using SC-LBM as being a parabolic surface,but was later challenged by Galindo-Torres et al.(2016) for two-phase displacement in multiple directions.Joekar-Niasar et al.(2010a,b)and Karadimitriou et al.(2014) numerically and experimentally validated its uniqueness for static cases.Still,both concluded the nonuniqueness of dynamic Pc-S-anwunder transient flow.For revisiting such a research interest,the dynamic Pc-S-anwsimulated in a multistep inflow/outflow test is shown in Fig.6a,with the projection on the Pc-S plane given in Fig.6b,and the projections of S-anwfor drainage and imbibition given in Fig.6c and d.

Fig.6.Dynamic Pc-S-anw: (a) Pc-S-anw in 3D coordinates,(b) The projection of Pc-S-anw in the Pc-S plane,(c) The projection of Pc-S-anw on S-anw for drainage,and (d) The projection of Pc-S-anw on S-anw for imbibition.

Fig.6a shows that a unique parabolic constitutive surface could not be fitted into all the data on Pc-S-anw.To confirm this,the three-dimensional(3D) parabolic surfaces were fitted into the Pc-S-anwfor one-to fifty-step inflow/outflow in Fig.7.Even though the fitting performance achieves high regression coefficients with 95%confidence bounds,their fitting parameters highly differ among different cases.If the one-and two-step inflow/outflow simulations are neglected,all other Pc-S-anwdata appear close and lay in a unique parabolic surface.The results in Fig.7 confirm the nonuniqueness of the constitutive surface of Pc-S-anwunder transient flow conditions (Joekar-Niasar and Hassanizadeh,2012a,b; Karadimitriou et al.,2014).When flow conditions reached a quasi-steady-state,the Pc-S-anwsurfaces finally converged into a unique constitutive relationship.

Fig.6b shows that only the one-and two-step inflow/outflow data deviated significantly from the steady-state Pc-S curves,while the rest of them deviated only slightly from the steady-state curves when pressure was varied on the boundaries in a step-wise manner.As shown by the curves of S-anwfor drainage and imbibition in Fig.6c and d,respectively,the S-anwcurves of multistep inflow/outflow in 5,10,20 and 50 steps all overlapped.By contrast,the S-anwcurves of multistep inflow/outflow in one-and two-step exhibited large anwunder highly transient inflow/outflow conditions.This further implies the stretching of the two-phase interface(meniscus) for fast drainage and the reverse of the meniscus for aggressive imbibition.Based on this,it appears that dynamic nonequilibrium effects under transient flow conditions can be alleviated if the pressure is increased in a smaller incremental interval with many more steps than five.When the pressure at boundaries was varied slowly enough,all Pc-S-anwdata overlapped into a unique parabolic surface.This could also be confirmed by Fig.6c and d,which show the uniqueness of S-anwfor both multistep drainage and imbibition with a number of steps of pressurization over five.

Note that the above simulations encountered an initializationrelated issue.Fig.6c and d shows a slight rise in anwfrom the outset of drainage and imbibition.This occurred because the parameters set in the fluid domain could not guarantee perfect equilibrium at the first time step.Although it took a quite short period to achieve the closest equilibrium condition,it was still inevitable to have a small abruption of anwnear the initial condition.

3.5.Flow regimes under transient multistep inflow/outflow

The dynamic Pc-S curves are usually related to either the first order of the temporal derivative of saturation (dS/dt)(Hassanizadeh et al.,2002) or the seepage velocity (q)(Wildenschild et al.,2001).As the dynamic contact angle of the meniscus is a function of the capillary number (Ca) (Sheng and Zhou,1992),the flow regime of transient multiphase flow yielding the dynamic Pc-S was also analyzed with respect to Ca.According to Eq.(20),Ca quantifies the ratio of the viscous effect to the surface tension.However,merely considering these two effects is not sufficient.Inertia needs to be considered as well.The Weber number (We=Re.Ca) quantifies the ratio of inertia to surface tension but neglects the viscous effect.Therefore,the Reynolds number (Re),which is rarely considered for the analysis of multiphase flow regimes,was used to estimate the inertia over the viscous effect using Eq.(19).The dynamic Pc-S curves of drainage and imbibition were simulated using multistep pressurizing boundary conditions,and are separately plotted against Ca and Re in Fig.8a and d.To better visualize the curves of S-Ca and S-Re,the Pc-S-Ca and Pc-S-Re curves are separately projected into the SCa plane in Fig.8b for drainage and Fig.8c for imbibition,and into the S-Re plane in Fig.8e for drainage and Fig.8f for imbibition.

Fig.7.Fitting the parabolic surface into Pc-S-anw:(a)One-step inflow/outflow,(b)Two-step inflow/outflow,(c)Five-step inflow/outflow,(d)Ten-step inflow/outflow,(e)Twentystep inflow/outflow,and (f) Fifty-step inflow/outflow.

Fig.8.Dynamic Pc-S against Ca and Re: (a) Pc-S-Ca in three dimensions,(b) The projection of the Pc-S-Ca for drainage into the S-Ca plane,(c) The projection of Pc-S-Ca for imbibition into the S-Ca plane,(d)Pc-S-Re in three dimensions,(e)The projection of Pc-S-Re for drainage into the S-Re plane,and(f)The projection of the Pc-S-Re for imbibition into the S-Re plane.

Fig.8b,e,c and f shows that since the commencement of inflow/outflow,both the Ca and Re,as functions of the seepage velocity(q),increased and finally decreased to zero at final equilibrium for oneand two-step inflows/outflows.Without calculating Re,it is impossible to determine whether inertia dominated the flow regime for fast drainage and imbibition.For the multistep outflow in 5-50 steps,Ca and Re fluctuated corresponding to step-wise increments in the fluid pressure along the boundaries.A comparison between inflow/outflow in a few(1-2)and many(5-50)steps of pressurization shows that when dynamic Pc-S further deviated from the quasi-steady-state case (in 50 steps),a highly transient multiphase flow occurred that was dominated by inertia.For these conditions,Re reached a value of 12,indicating a transition between the Darcy and Forchheimer flow regimes.In this case,the inertial effect was more significant than the viscous effect and surface tension.On the contrary,the inertia could be reduced when the number of steps of pressurization was increased at a smaller incremental interval.For instance,Re was less than five when the pressure was increased in five steps along the boundaries.In this case,Ca was less than 0.01.A smaller q thus ensured that the flow regime was dominated by the viscous effect or even surface tension when q was very small (i.e.Ca ＜0.001 and Re ＜1).

Due to the initialization issue mentioned above,as for the pressuring step in 5-50 steps,there are abruptions of q after the initial time frame,thus corresponding Ca and Re also present a sudden increase.A similar sudden increase in anwcan be observed in Fig.6c and d.Such an imperfect equilibrium for the initial condition needs to be solved in future work.

The results in Fig.8 show that the dynamic Pc-S curves generated by the multistep transient inflow/outflow simulation occurred because of the high pressurization of the fluid along the boundaries within a few steps (one to two).The inertia became dominant in such a highly transient multiphase inflow/outflow.When the boundary condition for pressure was adjusted sufficiently slowly(in 5-50 steps),the Pc-S curves under transient flow conditions approached the steady-state Pc-S curves.Moreover,the correspondingly lower Ca ＜0.05 and Re ＜2 indicate that the flow regimes returned to the state where either the viscous force or surface tension was dominant.Furthermore,as Ca and Re were functions of q,they reflected the unsteady acceleration q(t)missed in the momentum balance when formulating Darcy-type seepage equations for two-phase fluids.Therefore,the inertia,including unsteady and convective acceleration,should also be considered in the dynamic nonequilibrium seepage equations when inertia predominating the seepage flow under extreme boundary conditions(Liakopoulos,1964; Sposito,1980; Yan et al.,2021b; Chen et al.,2022; Yan et al.,2022c).

4.Conclusions

This study implemented the multiphase multicomponent SCLBM to investigate multistep and cyclic two-phase displacements in 2D porous media.The authors applied multistep and cyclic pressurizations on the boundaries of the simulation to replicate transient and steady-state multiphase inflows/outflows.Furthermore,the simulation setup enabled a pressure cell experiment as well as the point-wise measurement of capillary pressure(Pc) and saturation (S) on the profile of a porous media by applying multistep and cyclically vibrating hydraulic Dirichlet boundaries.The following conclusions were obtained:

(1) The results showed that the dynamic Pc-S could be simulated by multistep pressurization boundaries with large increments in pressure within a few steps.If the number of steps was increased with a smaller increment in pressure in each step,the Pc-S curves tended to be closer to curves of the steady-state case,with some deviations corresponding to the step-wise variation in the pressure along the boundaries.

(2) The fluid patterns,meniscus dynamics,phase trapping and local density around the two-phase interfaces between inflow/outflow simulations of cases involving a few steps(fewer than five) and more steps (more than five) were not precisely the same.This might be the reason why SC-LBM could simulate the dynamic Pcfor dynamic Pc-S curves under transient two-phase flow conditions.

(3) The cyclic inflow/outflow was generated using a triangularwave pressurizing boundary condition,such that the dynamic hysteretical Pc-S scanning loops could be simulated by SC-LBM.The dynamic hysteretical scanning curves also overshoot and undershoot out of the steady-state drainage and imbibition Pc-S curves,respectively.To the best of the authors’ knowledge,this is the first simulation of dynamic scanning curves using SC-LBM.

(4) The constitutive surface Pc-S-anwwas investigated based on the results of the multistep inflow/outflow simulations.The results showed that it would be impossible to determine a unique constitutive surface if the highly dynamic Pc-S curves under transient flow conditions were considered.However,with increasing the number of steps of pressurization with smaller increments in pressure,the differences between dynamic and steady-state Pc-S curves became smaller.For those cases,there were insignificant enlargements of the interfacial area during drainage and imbibition,thus the concept of the constitutive surface seems only suitable for such quasi-steady-state flow conditions.

(5) The dynamic Pc-S curves were investigated in terms of capillary number(Ca)and Reynolds number(Re).The results showed that significantly dynamic Pc-S curves were obtained in the transitional flow regimes between Darcy and Forchheimer,with Re ＞10 indicating that inertia dominated the flow.On the contrary,when Re ＜5,the dynamic Pc-S curves deviated only from the steady-state curves,with several fluctuations corresponding to multistep pressurization along the boundaries.

In principle,it is possible to conclude that the dynamic Pc-S curves are dominated by the boundary conditions of dynamic pressure from the perspective of a single REV scale.Significant overshoot and undershoot of capillary pressure for drainage and imbibition,respectively,occurred by varying the pressure on the top and bottom boundaries when using SC-LBM.Moreover,the greater the dynamic effects obtained under highly dynamic boundary conditions of pressure were,the higher the seepage velocities(q)were,which led to Ca ＞0.05 and Re ＞2.An analysis of the flow regimes indicated that inertia dominated flow over the viscous force and surface tension,and this reflected the unsteady acceleration q(t) missed in the conventional derivation of the momentum balance in the equations of two-phase seepage.Therefore,the theory of dynamic two-phase flow in porous media under transient flow conditions should be further developed by including not only the dynamic nonequilibrium Pc,but also the dynamic boundary conditions and inertial terms in the momentum balance in future work.

Finally,it is worth reflecting on the limitations of this study in comparison with other research on the multiphase LBM.First,SRT with unit density and viscosity ratios was used to reduce computational cost and enhance numerical stability.While many studies have used this setting in their SC-LBM,such a scenario is not universal in applications.Therefore,MRTs that permit high density and viscosity ratios will be implemented by the authors in future work to explore their dynamic effects using LBM.Second,the contact angle was set to zero with the Shan-Chen pseudo-potential to simplify the simulation setup.However,this is also not the case in most instances of natural two-phase displacement in porous media.This study focused on the effects of various pressure boundary conditions on the dynamics of capillary pressure and transient twophase flow in porous media.It would be informative to incorporate more effects,e.g.changes in wettability,higher density and viscosity ratios,and a variety of pore structures.Hence,the authors plan to extend multiphase schemes in Mechsys from the Shan-Chen formulation.Finally,this study applied only D2Q9 in quasi-2D due to constraints on computational power.The D3Q15 can be executed in three dimensions in future work as computational power improves.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The first author thanks the University of Queensland International Scholarship (UQI) for its support (Grant No.42719692).The simulations were implemented using the Mechsys open-source library on the Goliath high-performance supercluster at the University of Queensland,Australia.