Three‐dimensional Lattice Boltzmann simulation of gas‐water transport in tight sandstone porous media: Influence of microscopic surface forces

A three‐dimensional two‐phase Lattice Boltzmann model was developed and validated for the fundamental phenomena of gas‐water transport in tight porous media considering microscopic forces, namely, the electrostatic and solid‐liquid intermolecular forces. The thickness of the water film adhering to the porous media surface was calculated based on the principle of thermodynamic equilibrium and the gas‐water electrostatic potential energy. The Lattice Boltzmann model simulations focused on the effects of the microscopic forces and water film thickness on the gas‐water transport in tight porous media. The fluid flow in the porous media was simulated under different conditions for the characteristic length, interfacial tension, and displacement pressure gradient. The results confirmed that the gas‐water transport is strongly affected by the microscopic forces in tight porous media and that the transport process is accompanied by a high seepage resistance. The seepage resistance is more pronounced in the case of smaller pores because the adhering water film is thicker and more stable. During the transport process, the water film may disintegrate under certain conditions, which would enhance the gas flow in the porous media. However, it is difficult to effectively improve the gas flow by increasing the displacement pressure gradient under microscopic forces; this can be attributed to the accumulation of water or increased seepage resistance dictated by the interfacial tension, which occurs more frequently with smaller pores. This paper provides a promising new perspective for efficiently improving the Lattice Boltzmann model for transport trends considering microscopic forces.


| INTRODUCTION
The depletion of conventional gas and oil reserves has increased the attractiveness and strategic importance of unconventional energy resources such as shale oil and gas, tight sandstone oil and gas, and coalbed methane worldwide. Tight sandstone gas and shale gas reservoirs comprise ultratight porous media with pore sizes at the micrometer to nanometer scale, and the intrinsic matrix permeability and porosity are extremely low. China has abundant unconventional natural gas reserves; among these, the development of the tight sandstone gas reservoirs has the most practical significance because of their large size. 1,2 Zou et al 3 indicated that tight sandstone gas reservoirs in China have a pore throat radius of 25-700 nm and intrinsic matrix permeability of less than 0.1 × 10 −3 µm 2 under overburden pressure conditions. Such tight sandstone consists mainly of nanopores, which are the main passageway for fluid flow. 4 The most common mining technique is depletion; however, tight sandstone gas reservoirs always have high water saturation, 5 so the resistance to gas transport is considerably high, which reduces the gas recovery. Moreover, tight porous media have a complex microstructure; the microchannels and nanochannels reduce intermolecular collisions and enhance molecule-rock surface collisions, 6 and the surface of clay materials is usually negatively charged. Consequently, the surface effect is significant for the gas transport mechanism, particularly at a high Knudsen number (K n ). 7,8 Among all, the factors that influence the gas transport mechanism in tight porous media, the solid-liquid particle attachment is key. A reliable microscopic mechanical analysis of this factor is important for ensuring successful gas transport. In this study, the Lattice Boltzmann method (LBM) was used to develop a gas-water two-phase transport model that considers the influence of microscopic forces on the interface between the rock surface and fluid.
In the past decades, research on microscopic effects has focused on the gas flow characteristics in microchannels by investigating the drag coefficient of gas in micropipes. [9][10][11] Harley et al 12 studied the flow of gases in microtubules and found that forces had a considerably different effect on the flow at the microscale compared with at the macroscale. In general, the forces acting on a fluid are mass and surface forces. At the microscale, the effects of the surface forces are enhanced relative to the mass forces. In addition, the effect is more pronounced for regions closer to the rock pore surface. For example, if the length is reduced from 1 mm to 1 µm, the area and volume are, respectively, reduced to 10 −6 and 10 −9 times their original values. Moreover, surface forces such as the friction and viscous and interfacial tension forces increase 1000 times compared to forces proportional to the mass, such as inertia. Therefore, surface forces have a significantly increased effect on the flow characteristics at the micro-and nanoscales and cannot be ignored. Liu et al 13 analyzed the microscopic flow interface phenomena and flow boundary conditions and summarized the effects of the electrostatic force, van der Waals force, and spatial position force on the flow law of a fluid. The molecular forces were noted to play a key role in many related phenomena such as physical adsorption, surface infiltration, and film properties. Wang and Wu 14 studied the effect of an electric double layer caused by the electrostatic force, established an equation for fluid seepage in a microchannel under the action of the electric double layer, and analyzed the influence of the electric double layer on the microscopic flow law. Based on the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory, Derjaguin and Churaev 15 analyzed the relationship between the thickness of the adsorbed fluid film and pressure difference across the interface film considering the influence of the van der Waals force and diffusion double layer. Tuller and Dudley 16 and Tokunaga et al 17 applied the DLVO theory to studying the fluid film thickness in underground porous media connected with the atmosphere but did not consider that the film was trapped by the throat. In experiments, Stipp et al 18 observed that the water film adsorbed onto the surface of rock particles, and the adsorption area gradually increased. Numerical simulations based on molecular dynamics indicated that the adsorption of a monolayer liquid on the solid wall surface of a nanoporous channel caused the concave liquid surface to deviate from the ideal circular cross-section, 19 which was due to the reduced separation pressure from the interaction force between the wall of the microchannel and the fluid molecules. Li et al 20 analyzed the fluid saturation of a water-bearing shale gas reservoir, and its influence on the flow characteristics based on the thermodynamic equilibrium theory for the water phase. They established a model for calculating the thickness of the bound water film under the ultralow saturation condition. Their results showed that, during the formation of a water-bearing shale gas reservoir, the vaporization and liquid-carrying phenomena were the main reasons for the ultralow water saturation. In addition, some experimental results indicated that the saturated water core was displaced by the gas phase to an irreducible water saturation state (23.90%-25.43%), and the water saturation was further reduced to 12.4%-14.75% by dry gas flooding. 21,22 Therefore, accurately characterizing the thickness of the stable water film adhering on a sandstone wall is a key issue to developing a fluid flow law for micro-and nanoscale porous media. The effect of microscopic surface forces between the water phase and solid rock wall can be reflected by the thickness of the stable adhering water film. In general, negative charges are presented on the sandstone wall, and a diffusive double layer forms at the solid-liquid interface after they contact the formation water. A strong molecular force exists between the surface of the rock particles and the water phase, and a stable adhering water film often forms on the wall surface.
Another important aspect of determining the microscopic gas-water flow is the research methodology. Depending on the Knudsen number, the gas transport in porous media consists of four flow patterns: viscous flow, slippage flow, transitional flow, and free molecular flow. 23,24 For conventional reservoirs with large pore and throat sizes, the gas transport mechanism corresponds to viscous flow, which can be characterized by Darcy's law. For tight porous media, the Knudsen diffusion cannot be ignored, and Darcy's law is not valid [25][26][27] ; consequently, conventional fluid flow simulation methods are no longer applicable. Several scholars 28,29 have performed experiments considering natural or artificially bonded porous media and demonstrated that the particle size and porosity have a significant effect on the fluid flow. However, experimental methods are rather time-consuming and expensive. Furthermore, the experimental conditions for fluid flow in ultratight porous media are extremely challenging to achieve. Thus, numerical simulations based on solving the Navier-Stokes (N-S) equation are desirable for studying problems involving multiple mobility mechanisms at the micro-and nanoscales. The LBM is a comprehensive approach to solving the N-S equation. 30 The simulated object does not correspond to a large number of individual fluid molecules but instead consists of fluid particles (ie, molecular microclusters) that are sufficiently large at the microscale and sufficiently small at the macroscale, which considerably reduces the number of particles. 31,32 The LBM is suitable for fluid system modeling and simulation at the mesoscale and has seen extensive developments worldwide in recent decades. The LBM falls between the macrocontinuum hypothesis and microscopic molecular dynamics simulation, 33 so it combines the advantages of multiscale fluid flow simulations. As a starting point, the LBM directly discretizes the model based on the three conservation laws of matter: mass conservation, momentum conservation, and energy conservation. It also combines the theory of molecular motion and statistical theory to bridge the macroscale and microscale and the continuous and discrete aspects so that fluid motion can be characterized from a new perspective. As given in Table 1, the overall motion characteristics of particles are expressed by a distribution function, which is the key to the LBM and represents the collective behavior of a large number of particle groups. This study explored the gas and water flow characteristics in tight sandstone porous media under microscopic force conditions. Different types of fluid transport mechanisms that do not follow the continuum hypothesis were considered with the LBM. The general expressions of microscopic forces are at the macroscale; however, discrete forms at the mesoscale are required for the LBM to be applicable. To address this challenge, another thermodynamic method and Darwin's law were used to simulate the surface evolution of microscopic forces.
In recent years, rapid advances with the LBM have led to a series of methods that include multiphase fluid flows, high-speed flows, and so on. Hu et al 34 40 used the computational fluid dynamics method and the experiment method in order to analyze the particle deposition under various conditions and concluding that Boltzmann equation can well describe the relationship between gas Reynolds number and particle deposition in the rectifying plate system.
The aim of this study was to develop a three-dimensional two-phase Lattice Boltzmann model for exploring flow in tight porous media at the micro-and nanoscales with various pore and throat complexities while considering the influence of microscopic surface forces. Simulations of the gas-water transport confirmed the applicability of the gas-liquid interfacial tension and gas-liquid two-phase separation. Furthermore, the van der Waals force and electrostatic force from the rock particle wall were considered to quantify the effects of some relevant influencing factors on the breakthrough of the gas at the outlet. The mechanisms of gas-water transport in tight porous media were also clarified. The remainder of the paper is organized as follows: Section 2 briefly introduces the pseudopotential Lattice Boltzmann model of gas and water and describes the model validation for the gas-water interfacial tension and gas-water two-phase separation. Section 3 describes the calculation of the thickness of the stable water film adhering on the rock surface and the establishment of the two-phase Lattice Boltzmann model considering microscopic forces. Section 4 discusses the breakthrough time of the gas front with different influencing factors and the gas-water transport mechanism in detail. Finally, Section 5 discusses the results and conclusions.

| Lattice Boltzmann model
The Bhatnagar-Gross-Krook (BGK) evolution equation for a single relaxation time is an approximation of the continuum Boltzmann equation and can be discretized in the time, space, and velocity domains to obtain a fully discrete LBM equation. 41 The present three-dimensional two-phase LBM simulation is based on the BGK evolution equation and Shan-Chen model. 42,43 The fluid in the pore is considered to comprise two phases: gas and liquid. The gas phase is natural gas under the formation condition and contains multiple components, and the liquid phase is the formation water. In this study, the change in components during gas motion was not considered. For a two-phase gas-water flow in tight porous media, the probability density distribution functions of the discrete gas and water particles are given below: where f g is the probability density distribution function of the gas, f w is the probability density distribution function of the water, x is the space position vector, e i is the lattice velocity vector in the ith direction, and t is the time. τ g and τ w are the gas and water relaxation times, respectively, and can be determined by considering the fluid kinematic viscosity. f g(eq) and f w(eq) are the gas and water equilibrium distribution functions, respectively, and ∆t is the unit time of the LBM. In general, finding a numerical solution to the LBM equation involves two major steps: migration and collision. Discrete fluid particles can only migrate along the lines of the lattices and can only move into the neighboring lattices according to a fixed velocity scheme within each time step. During the collision step, the probability density distribution function relaxes toward its equilibrium state at a given relaxation time. The equilibrium distribution functions of the discrete gas and water particles are given below: where ρ g and ρ w are the mass densities of the gas and water components, respectively; u g(eq) and u w(eq) are the equilibrium velocities of the gas and water, respectively; c s is the lattice sound velocity; and w i is the weighting factor for the velocity distribution. Figure 1 shows the scheme of the D 3 Q 19 model, which contains 19 discrete velocities in the three-dimensional space. For these 19 directions, w i is 1/3 at i = 0, 1/18 at i = 1-6, and 1/36 at i = 7-18. Based on the definition of the Lattice Boltzmann model, the densities of the gas and water phases at the macroscale are given as According to the Shan-Chen pseudopotential model, the strength of the interaction between the two phases is characterized by the interaction potential energy between particles, which also governs the phase interface capture and tracking.
The surface force F x,x ′ acts on the σ particles from neighboring particles. To maintain consistency with the expression of the pseudopotential model, σ is used to represent the gas and liquid phases rather than the previous two equations. The surface force is expressed as follows: where ψ σ is the effective density of component σ, which acts as a density-dependent potential function or effective mass related to the density of particle σ. To make the numerical results of the model consistent with those of the Maxwell thermodynamic theory of the equation of state, the effective density can be taken as ψ σ = ρ 0 [1 − exp(−ρ/ρ 0 )]. ρ 0 is the normalized constant of the density and generally has a constant value of 1. Ĝ represents the interaction strength between the neighboring particles σ and, which can be the same or different types of particles. Ĝ = 0 is in a nonadjacent lattice, while Ĝ = G is in an adjacent lattice. G is the coefficient of interaction strength. x and x′ represent two adjacent positions. When x → x′, and V is the interaction potential function of gas and water.
This indicates that the external force term of the pseudopotential model is given by where F i denotes the gas-liquid and solid-liquid interaction surface forces (ie, external forces), and x + e i {i = 1, 2, …, N} corresponds to N lattice points centered on the x lattice. The strength of the interaction between different particles is represented by the interfacial tension between different types of particles. The extraforces are contained in the equilibrium state distribution functions, which are used to characterize the interactions between the gas and water. Therefore, the corresponding evolution equation of the distribution function is given by When no external force exists, the local momentum is conserved in the collision process, so no interaction occurs between the gas and water. Thus, the average lattice velocity is given by where u σ denotes the flow velocity of component σ.
The velocity of the equilibrium state of the particles under an external force is given by where u σ(eq) denotes the equilibrium velocity of component σ.
The fluid density can be calculated as where u is the average velocity of the mixed fluid before and after a collision, and = ∑ is the mixing density of the gas and water. The derivation process and more details concerning the pseudopotential Lattice Boltzmann model are given in previous studies. 44,45

| Gas-liquid interfacial tension
The interfacial tension between the gas and liquid causes the two phases to separate. The force between the gas and liquid (ie, surface tension) becomes an additional force in Equation (10) at the gas-liquid interface. The spurious velocity is an important parameter that restricts the accuracy of the LBM for multiphase fluids. The pseudopotential model does not directly provide the surface tension during the construction process, which influences the gas-liquid flow in such tight porous media. To avoid simulation errors when the pseudopotential model is used for a two-phase fluid separation simulation with a high density ratio, the gas-liquid interfacial tension was characterized by a modified form of the pressure tensor: This was proposed by Hu, 46 who focused on adjusting the pressure tensor of the pseudopotential Lattice Boltzmann model to reduce the influence of the spurious velocity. Here, P is the pressure tensor, A is the weight coefficient, κ is the interfacial tension coordination coefficient, and I is the unit tension perpendicular to the interface. The interface pressure perpendicular to the interface can eliminate κ by simplification, while the interfacial tension parallel to the interface can cancel out the terms related to I. To ensure the stability of the model, A is usually −0.5. The value of κ lies in [0, 0.99]. Hu 46 analyzed the spurious velocity of the improved pseudopotential model and showed that the spurious velocity decreases with increasing temperature. In the present study, the maximum spurious velocity in the simulations was only 1.56 × 10 −6 and thus had a little influence on the calculation results.
After simplification, the interfacial tension can be calculated as follows: where P it is the interfacial tension. For each gas-water interfacial tension, there exists a corresponding unique particle interaction strength; this can be used to characterize different interfacial tensions. Therefore, the relationship between the gas-liquid interfacial tension and interfacial tension coordination coefficient is as shown in Figure 2.
The force between the two phases at the gas-liquid interface is in the form of interfacial tension and is equal to the gradient of the gas-liquid interaction surface forces in Equation (10), which is given by ∇F i = P it . Thus, the interfacial tension at the gas-liquid interface can be defined as To improve the computational efficiency, an explicit method is utilized. When the probability density distribution gradient is calculated for a new time step, the gradient for the previous time step is used. The difference equation is expressed as follows: Generally, the value of the lattice space quantity is arbitrary. If the lattice scale is converted to the equivalent macroscale, the macroscopic Reynolds number for each gas and liquid phase Re σ = v σ d/γ σ needs to be equal to the dimensionless criterion number Re σ = u σ H/υ σ in the lattice flow field, where v σ and d σ denote the macroscopic characteristic velocity and length, respectively, of the σ phase. The corresponding equations for the conversion coefficients of the fluid parameters at the macroscale and lattice scale are presented in Table 2.
The Mach number is a measure of the fluid compressibility and is given by The velocity of the gas-liquid two-phase flow in a tight reservoir is slow, and its Mach number is less than 0.3 when the parameters are changed. Therefore, the two-phase flow can be considered as an incompressible fluid.

| Gas-liquid two-phase separation simulation
The two-phase pseudopotential Lattice Boltzmann model was developed in a computational program. Before it was used for relevant numerical studies of the gas-liquid transport in porous media, the model was validated for two test cases: gas-liquid interfacial tension and gas-liquid two-phase separation. According to Laplace's law, when a lattice system is in equilibrium, the two-phase static pressure difference ∆p between the contact surfaces is proportional to the interfacial tension σ between the fluids and inversely proportional to the interface radius of curvature r: In this case, the model had 300 × 300 × 300 lattices. The inlet and outlet were located on the left and right sides, respectively, of the model and subjected to the constant pressure boundary condition with a low pressure gradient of 0.01 MPa/ cm (7.33 × 10 −4 in the Lattice Boltzmann model) so that the gas-liquid separation could be observed. In addition, other boundaries (ie, the upper, lower, front, and back boundaries) corresponding to the directions perpendicular to the mainstream were set according to the periodic boundary hypothesis.
Because of the presence of the water phase, the boundary slippage effect could be ignored. 47 The kinematic viscosity coefficient υ was − 0.5 c 2 s Δt, and the relaxation time of different fluids can be inversely determined from the given fluid viscosity. Thus, the gas density was set as 1, and the water density was 1000 under the surface condition with a ratio of 1.37 under the underground condition (110 degrees Celsius, 30 MPa). The kinematic viscosities of the gas and water were 0.02 and 0.2, respectively. The gas-water interfacial tension was set to 0.003. All of the above parameters were dimensionless lattice units. To test Laplace's law, a series of LBM simulations was performed with different values for the interface radius of curvature r, as shown in Figure 3.
To simulate the gas-water separation, a density disturbance of 1% was first applied to the density of the gas-water simulation system. This was used to control the fluid density of each lattice to float within the preset range (error of 1%), which was equivalent to each lattice being filled with a large number of very small gas bubbles. Bubbles smaller than the lattice gradually merged into larger bubbles. The gas-water flow was simulated according to the two-phase flow lattice model described above. Figure 4 clearly shows the gas-water separation process at the middle section and the corresponding iterative time steps. The gas distribution position and content are indicated in terms of the ratio of the gas density in each lattice to the density when the lattice was filled with gas. The red region represents the pure gas phase, the blue region corresponds to the pure water phase, and the transition colors from red to blue represent different gas contents in the lattice. The evolution time steps were 100, 800, 2000, 5000, 8000, and 20 000 (ie, 7.01 × 10 −6 , 5.61 × 10 −5 , 1.40 × 10 −4 , 3.51 × 10 −4 , and 1.40 × 10 −3 seconds, respectively).
As shown in Figure 4, at 100 time steps (7.01 × 10 −6 seconds), the gas-water separation was in the initial stage with very small bubbles and no lattice for pure gas or liquid. At 800-8000 steps (5.61 × 10 −5 -3.51 × 10 −4 seconds), the gas and water gradually separated over time. Under the action of the interfacial tension, the gas gradually coalesced and formed bubbles, which eventually became larger. A small number of small bubbles did not fully integrate into large bubbles. At 20 000 steps (1.40 × 10 −3 seconds), almost all of the small bubbles merged into large bubbles, and these bubbles merged with each other. Only a few bubbles remained in the whole model. The pseudopotential model was able to separate the gas and liquid phases, but there was no solid phase. Tight porous materials have a water film at the solid-liquid interface and microscopic surface forces. Therefore, the proposed model requires further enhancement for practical applicability to tight porous materials.

WATER FILM AND MICROSCOPIC SURFACE FORCE IN MICROPORES
Tight sandstone has strong hydrophilicity and high irreducible water saturation in tight gas reservoirs. When the water and gas phases coexist in tight porous materials, a water film with a certain thickness and stable adsorption is usually presented on the walls of the porous media, which affects the mass transfer of the gas. 48 In micropores and nanopores, the throat structure for gas-water two-phase seepage in tight sandstone causes the solid-liquid microscopic surface forces between the solid rock wall and fluid water phase to be stronger than the gas-solid interfacial force for single-phase gas seepage. Therefore, the effects of the adhering water film and microscopic surface forces of the water phase must be considered in the Lattice Boltzmann model of the gas-water two-phase flow in a tight reservoir.

| Thickness of the stable adhering water film
The solid-liquid interfacial force between the walls of the tight porous material and water phase not only hinders the flow of the fluid water phase but also hinders the water relative to the gas flow. This interaction between the solid, liquid, and gas is reflected by the water film adsorbing stably on the rock wall. According to the maximum entropy principle, a thermodynamic equilibrium exists between the liquid water and water vapor that ensures the stability of the water film. Considering that the water phase is incompressible and the temperature of the tight sandstone porous media remains constant, the law of thermodynamics and Y-Laplace formula can be used to define the influence of the unit pressure difference on the molar free energy of the water phase 20 : where V w = M w /ρ w is the molar volume of water (cm 3 /mol), σ is the gas-water interfacial tension (N/m), and R 1 and R 2 are the radii of the main curvature. The pores of the porous media are assumed to be spherical, so the additional pressure difference at the gas-water interface can be expressed as ∆p = σ (1/R 1 + 1/R 2 ) = 2σ/r, which is the pressure difference between the water film and bubble (MPa). r is the radius of curvature of the inward water film surface; in this study, it was taken as the pore radius (m). Furthermore, the change in free energy of the water phase can also be expressed as the change in the vapor pressure. According to the G-M binary gas isothermal adsorption

Gas content (%)
model, 49 the adsorption potential energy to form 1 mol of water vapor from liquid water is given by where p 0 is the partial pressure of the water vapor on the bending water surface (MPa), p is the corresponding saturated vapor pressure on the gas-water surface (MPa), and the universal gas constant is R = 8.314 J/(mol K). According to the Kelvin formula, when the gas and water phases are in stable state in a porous media (ie, ∆G wl = G wg ), then Consequently, the pressure difference at the gas-liquid interface is Accordingly, the relative humidity (RH) of the natural gas is given by the following relation: where the relative vapor pressure p/p 0 is related to the pore radius r.
In tight porous sandstone reservoirs, because the nanopores are the main space for fluid storage and transport, the effects of microscopic forces on the fluid transport behavior should be considered. For water particles, the water film tends to adhere stably to the solid wall because of the interactions between the microscopic surface forces of the solid wall, water film, and gas. The pressure generated on the water film during this process is called the separation pressure Π 15 which is a function of the water film thickness h. The separation pressure usually consists of the van der Waals force Π m , electrostatic force Π el , and tectonic force Π s ; each is also a function of the water film thickness h. In addition, Π = Π m + Π el + Π s .
Because the main component of sandstone is quartz, which has strong hydrophilic properties, Masayuki's argument 20 was used to consider the van der Waals force and electrostatic force for the analysis on the effects of the microscopic forces on fluid transport behavior in tight porous media. The action of the diffused electric double layer generates the electrostatic force Π el (h); the double-layer model established by Langmuir 50 can be used to describe the water film thickness adhering on the surface of the media at a high electric potential energy: where ε r is the relative permittivity of water (F/m), ε 0 is the vacuum permittivity (F/m), k B is the Boltzmann constant (J·K −1 ), e is the absolute value of the electron charge (C), and Z i is the valence state.
The van der Waals force is generated by the dipole interaction between the solid wall, water film, and gas and can be expressed as 51,52 where A s-l-g is the Hamaker constant of interaction between the solid wall and gas through the intermediate water film and is generally approximated as 10 −21 J (7 × 10 −21 J in this study).
According to the thermodynamic equilibrium principle and G-M binary gas isothermal adsorption model, 28 when the water film and gas reach an equilibrium, the separation pressure is related to the relative vapor pressure as follows: where V w is the molar volume of water, the universal gas constant is R = 8.314 J/(mol K), and T is the Kelvin temperature.
Because the tectonic force Π s can be neglected at the scale of tight sandstone pores, Equations (24)- (26) can be used to obtain the following relationship: Equation (27) is used to calculate the thickness of the stable adhering water film considering the combined action of the van der Waals force and electrostatic force between the formation water and rock wall. According to related data from Mattia et al, 53 under stratigraphic conditions (110°C, 383.15 K), the water permittivity ε r is 80 F/m, vacuum permittivity ε 0 is 8.85 × 10 −12 F/m, the molar volume of water V w is 18 × 10 −6 m 3 /mol, and the gas-water surface tension σ is 33.86 × 10 −3 N/m under formation temperature conditions. Figure 5 shows the calculation results for the stable adhering water film thickness for different pore throat radii. The thickness of the stable adhering water film clearly increases with the pore radius. When p/p 0 approaches the critical point for stable coexistence between gas and water, the water film thickness increases rapidly, which forms capillary water in situ. When the relative humidity of the gas is fixed, a smaller pore throat corresponds to a greater thickness of the stable adhering water film, greater seepage resistance, and greater difficulty with gas migration. These results further explain why gas cannot be effectively displaced and is trapped in the corners of the porous media and tiny capillaries.

Boltzmann model considering microscopic surface forces
Because of the presence of the water film, the model only has a solid-liquid interface and no gas-liquid interface. The two-phase Lattice Boltzmann model considers the van der Waals force and electrostatic force between a solid surface and the water film adhering onto it. When the Lattice Boltzmann model is set up, as the introduced microscopic forces act on the surface of the porous media, the distribution function and evolution equation of the fluid (ie, water phase) at the boundary are corrected, and the additional external forces in the model are discretized. Based on Equations (23) and (24), the microscopic forces can be expressed as follows: By discretizing the microscopic forces at the solid-liquid contact layer, the corresponding equation is given below: where the additional external force term F i is a function of F, and f eq i is a discrete equilibrium distribution function that depends on the density ρ w and mass velocity (ρu) w . The pore size can be used to obtain the thickness of the water film adhering onto the surface. The Chapman-Enskog multiscale expansion can be applied to transform microscopic equations into macroscopic equations: where a is the acceleration caused by the van der Waals force and electrostatic force on the solid-liquid interface.
The velocity distribution of the water phase can be expressed as follows: The model described above is the two-phase Lattice Boltzmann model for gas-water transport in tight porous media under the effects of microscopic forces established in this study. The model is based on an analysis of the mechanics for the stable adherence of a water film on the media surface, which can be used to describe the two-phase gas-water transport behavior in microscopic tight porous media.

| RESULTS AND DISCUSSION
Numerical analysis was performed with the established Lattice Boltzmann model of the fluid transport considering microscopic forces. The characteristics of the water film thickness on the tight porous surface and influence of the microscopic forces on the two-phase gas-water transport behavior in tight porous materials were investigated.

| Lattice Boltzmann model stability test
To validate the model, stable conditions in the real world were assumed. The computational domain was set to 100 × 100 × 100 lattice nodes in the x, y, and z directions. The six boundaries of the model comprised the solid phase, the center of the model comprised a gas bubble with a diameter of 30 lattice nodes, and the rest was filled with water. The temperature of the lattice system was 110°C, and the pressure was 30 MPa. The gas kinematic viscosity was 2.60 × 10 −7 m 2 /s (0.2 in the model), while the water kinematic viscosity was 2.60 × 10 −8 m 2 /s (0.02 in the model). The density distribution and spurious velocity distribution of the model section are shown in Figure 6. The gas-liquid density ratio was 1.37 under the formation condition. The spurious velocity in the test was concentrated near the gasliquid interface, and the maximum spurious velocity was only 1.50 × 10 −6 . In contrast, the velocity in the model was generally greater than 0.01.

| Analysis of the influencing factors
Because the flow of the water phase (wetting phase) is affected by the van der Waals force and electrostatic force from the rock particle surface, the fluid flow law is influenced by the surface tension between the gas and water phases and the displacement pressure gradient. In addition, because the pore throat is at the micro-and nanoscales, the scale effect is more prominent. In this study, a series of LBM simulations was performed to analyze the influencing factors for the gas-water transport. The developed LBM gas-water two-phase flow model was used to simulate the process of gas flooding the water. The computational domain was set to 1000 × 100 × 100 lattice nodes in the x, y, and z directions, respectively. The simulation was performed with a regular circular channel. The voids in the y-z section were arranged in a circular shape surrounded by solid boundaries, as shown in Figure 7.
The pseudopotential model was used to calculate the gas-liquid interaction force, and the additional force was calculated with Equation (17). The additional external force from Equation (29) was used to set the solid-liquid boundary conditions. All of the lattices were set to be filled with water; the temperature of the lattice system was 110°C, and the pressure was 30 MPa. The gas kinematic viscosity was 2.60 × 10 −7 m 2 /s (0.2 in the model), while the water kinematic viscosity was 2.60 × 10 −8 m 2 /s (0.02 in the model). Given a certain displacement pressure, the gas was injected from the inlet end of the model, and the breakthrough of the gas at the outlet of the model was taken as the condition for the simulation to terminate.

| Characteristic length
The three-dimensional Lattice Boltzmann model was used to simulate the two-phase gas-water flow with four characteristic lengths (20,50,100, and 200 nm) as well as with and without the electrostatic force and van der Waals force at the wall boundary. The displacement pressure gradient was set to 0.  characteristic length could be used to calculate the thickness of the stable adhering water film and obtain the action potential of the solid-liquid wall surface. Based on the dimensional analysis and scale conversion, the time step at the grid scale could be converted into the leading-edge breakthrough time in seconds at the macroscale. For example, at a characteristic length of 200 nm, the similarity criterion indicates that the length in the x direction of the lattice model can be converted to 2000 nm at the macroscale. Consequently, the breakthrough time of the gas leading edge as the gas flooded the water was the time required for the gas leading edge to pass through the plate with a length of 2000 nm. The method for calculating the distance covered by the gas leading edge before its breakthrough is described below.
The simulation results are shown in Figure 8.
The results indicated that a smaller characteristic scale corresponds to a lower gas flow rate, and longer time required for the leading edge to break through at the model outlet during gas flooding. When the interaction effect of microscopic forces on the water phase and rock wall was considered, the gas transportation became more difficult, and the flow velocity was lower. Therefore, it took longer for the gas leading edge to break through the small pore throat. As the characteristic scale increased, the influence of the microscopic forces on the gas transport capacity gradually decreased. At a characteristic scale of 200 nm, their influence on the gas transport capacity was sufficiently small that it could be ignored. The microscopic forces between the rock wall and water phase acted on the water body to form a stable adhering water film with a certain thickness. It has previously been proven that a smaller scale corresponds to a more stable adhering water film and greater thickness. A tight sandstone gas reservoir has small pore sizes, and the stable adhering water film on the rock wall has a greater retention effect on the gas transport. Therefore, the gas transport is subjected to greater resistance at extremely small characteristic scales ( Figure 9).

| Interfacial tension
The interfacial tension between gas and water is the difference in free energy at their interface, which directly reflects the ability of the water phase to hinder the gas flow. The correlation in Figure 2 shows the interfacial tension can be converted to the interfacial adjustment parameter to change the other model parameters in the Lattice Boltzmann model. The characteristic scale was set to 50 nm. The displacement pressure gradient at the inlet and outlet was set to 0. gas leading edge. This decreased the sweep efficiency and increased the irreducible water saturation. When the solid-liquid microscopic forces were considered, the interfacial tension was higher, and the gas flow was more hindered. The microscopic forces between the solid rock wall and water phase indirectly increased the difficulty of gas transport through the stable adhering water film. Regardless of the effect of any of the microscopic forces, the water film was the key to obstructing the flow of gas in the tight porous media. To improve the gas transport in micropores and nanopores, the restraint of the water film needs to be broken to achieve efficient development of tight gas reservoirs; the most effective means of doing so are fracturing.
The characteristic scale was set as 50 nm, and the gas-water interfacial tension was set to 0.04 N/m (5.90 × 10 −4 in the model. As the gas displaced the water, a higher displacement pressure gradient corresponded to the gas leading edge breaking through earlier at the outlet of the model. When the solid-liquid microscopic force was considered, the gas leading-edge transport was delayed at the same displacement pressure gradient, and it took a long time for the leading edge to break through. Theoretically, a greater displacement pressure gradient corresponds to a higher fluid transport velocity. However, Figure 10 shows that increasing the displacement pressure gradient did not significantly improve the gas transport in the water-bearing porous media. Under the action of the electrostatic force and van der Waals force between the water phase and solid rock wall, it was more difficult to drive the water with the gas. This was because increasing the displacement pressure gradient aggregated part of the water phase that was originally the stable adhering water film, and a larger amount of the water phase participated in the flow to generate additional resistance to the gas flow. Therefore, when exploiting tight sandstone gas reservoirs, the appropriate displacement pressure gradient should be selected based on the geological characteristics of the target field. Indiscriminately increasing the displacement pressure difference is difficult and costly and does not necessarily lead to the expected development.

| Gas-water transport mechanism in tight porous media
Numerical analyses were performed on the gas-water twophase transport behavior in tight porous media. Figure 11 shows scanning electron microscope (SEM) images of a tight sandstone core obtained from the Sulige Gas Field in China. The image resolution is 1.9 nm. The porosity of this tight sample was 8.62%, and the overburden permeability was 0.0937 × 10 −3 µm 2 . Figure 12 shows the corresponding microstructures of the tight porous media sample obtained by three-dimensional reconstruction. In this case, the two-phase Lattice Boltzmann model of the gas-liquid flow under microscopic forces was used to simulate the process of the gas flooding the water. The main focus was on the transport behavior of the fluid in the porous media at the micro-and nanoscales. The computational domain was set to 300 × 300 × 300 lattice nodes in the x, y, and z directions, respectively. Similar to the analysis on the influencing factors, the pseudopotential model was used to calculate the gas-liquid interaction force. The additional force was calculated with Equation (17), and the additional external force from Equation (29) was used to set the solid-liquid boundary condition. The temperature of the lattice system was set to 110°C, and the pressure was set to 30 MPa. All lattices were set as filled with water. The displacement pressure gradient was set to 0.6 MPa/ cm (1.76 × 10 −4 in the model), while the pressure was constant during the entire displacement process. The process of the gas displacing the water is triggered by disturbing the density at the inlet end. The gas kinematic viscosity was 2.60 × 10 −7 m 2 /s (0.2 in the model), while the water kinematic viscosity was 2.60 × 10 −8 m 2 /s (0.02 in the model). The gas-water interfacial tension was set to 0.03386 N/m (0.01 in the model). Figure 13 shows the fluid distribution in the initial state (water saturation of 100%). Figure 14 shows the gas-water distribution after the end of the displacement process, (ie, stable distribution). The gray part represents the solid wall, the blue part represents the gas phase, and the yellow part represents the water phase.
The simulation results showed that the gas front achieved a breakthrough at 89 800 time steps (3.15 seconds). The transport of gas and water reached a stable state at 120 000 time steps (4.21 seconds). These results indicate that some optimal channels exist for gas seepage in tight porous media under water-bearing conditions. The porous media are so tight that, once the optimal transport channels have formed, all other connected pores are auxiliary paths. Consequently, the gas and water transport quickly reaches a stable state after the breakthrough of the gas front. Figures 15 and 16 show the gas and water distributions for different model sections at different positions. Figure  15A shows the gas and water distribution at the position of Z = 240 along the Z direction when the gas front broke through at the exit of the model. Circle A in Figure 15A indicates that the water was preferentially displaced in the larger pores. When the gas leading edge broke through, the gas phase was mainly at the center of the larger pores, and the water phase was on the solid pore surface in the form of a thin water film. As shown in circle B, the thickness of the water film accounted for a larger proportion of the pore space for smaller pores. This is consistent with the results obtained from the calculation model for the stable adhering water film as described above. With the displacement, the distributions of the different fluids tends to become stable. Circle C in Figure 15A,B shows that the water can be driven into some smaller pores. However, because of the stronger microscopic forces acting on the solid surface and water, a thicker adhering water film was on the solid surface, which significantly reduced the fluid transport capacity in the porous media. Furthermore, when the gas front broke through, part of the water in the larger pores was still not driven because of local retention of the gas. Thus, a continuous transport channel did not yet form, such as circle D in Figure 15A. Figure 16 illustrates the stable state of the gas and water distribution at the position of Z = 110 along the Z direction. The irreducible water in the corners and tiny pores could not be effectively driven throughout the model. Overall, the fluid transport capacity was lower in the tight porous media. Because of the dual effects of the gas-water surface forces and solid-liquid microscopic forces, the water phase stably adhered to the solid rock surface, which reduced the cross-section of the gas transport channel, while the two-phase gas-water transport in the tight porous media increased the seepage resistance. The gas was transported at a lower speed in smaller pores because it tended to be trapped by the thicker water film. Therefore, the transport capacity of gas is considerably reduced in aqueous tight porous media.

| CONCLUSIONS
In this study, the pseudopotential LBM was coupled with the electrostatic force and solid-liquid intermolecular force to investigate the two-phase transport mechanism of gas and water in tight porous media. The influence of microscopic forces was used to calculate the thickness of the stable water film adhering on a rock surface, and the effects of varying values for the characteristic length, interfacial tension, and displacement pressure gradient on the gas-water transport mechanism were examined. The microscopic forces were found to considerably influence the water film on the rock surface and the trends for the gaswater flow. The main conclusions of this study are as follows: 1. The thickness and stability of the water film adhering on to the rock wall surface can be attributed to the interactions of the solid-liquid molecular force and electrostatic force. A smaller pore throat corresponded to a greater thickness of the water film and a higher transport resistance. This aspect was why the gas in the corners and very small capillary tubes could not be displaced and was retained in the pore channel. In contrast, as the pore throat size increased, the thickness of the water film decreased.

2.
The simulation results showed that the presence of microscopic forces increased the difficulty of fluid transport in the region adjacent to the wall surface and affected the gas and water flow by thickening the stable adhering water film and directly impeding the fluid transport. When the influence of the microscopic forces on the actual rock wall surface is considered, the fluid transport capacity is considerably different from that without the microscopic forces, particularly around the region adjacent to the rock surface. 3. The analysis on the influencing factors for the breakthrough time of the gas front indicated that the effects of the microscopic surface forces could not be ignored when the characteristic scale was less than 200 nm. A smaller characteristic scale corresponds to a greater seepage resistance. Greater gas-water interfacial tension corresponds to greater bound water saturation and more difficult gas transport. Increasing the displacement pressure gradient does not significantly improve the gas transport capacity in water-bearing tight sandstone gas reservoirs. 4. Under the combined actions of the gas-water surface force and solid-liquid microscopic forces, the water phase stably adheres to the rock wall surfaces of the pore and throat, which reduces the seepage section area. The movable water is driven out under the action of the differential pressure, which forms a two-phase gas-liquid flow and enhances the seepage resistance. The gas transport velocity is lower in smaller pores with a thick film of bound water, where the gas is easily trapped by the water phase. These factors comprehensively reduce the transport capacity of the gas phase in water-bearing tight sandstone gas reservoirs.