# Effect of Hydrate on Gas/Water Relative Permeability of Hydrate-Bearing Sediments: Pore-Scale Microsimulation by the Lattice Boltzmann Method

Authors: Xin Xin; Bo Yang; Tianfu Xu; Yingli Xia; Si Li
Journal: Geofluids
 (2021)

---

## Abstract

As a clean energy source with ample reserves, natural gas hydrate is studied extensively. However, the existing hydrate production from hydrate deposits faces many challenges, especially the uncertain mechanism of complex multiphase seepage in the sediments. The relative permeability of hydrate-bearing sediments is key to evaluating gas and water production. To study such permeability, a set of pore-scale microsimulations were carried out using the Lattice Boltzmann Method. To account for the differences between hydrate saturation and hydrate pore habit, we performed a gas-water multiphase flow simulation that combines the fluids’ fundamental properties (density ratio, viscosity ratio, and wettability). Results show that the Lattice Boltzmann Method simulation is valid compared to the pore network simulation and analysis models. In gas and water multiphase flow systems, the viscous coupling effect permits water molecules to block gas flow severely due to viscosity differences. In hydrate-bearing sediments, as hydrate saturation increases, the water saturationSw between the continuous and discontinuous gas phase decreases from 0.45 to 0.30 while hydrate saturation increases from 0.2 to 0.6. Besides, the residual water and gas increased, and the capillary pressure increased. Moreover, the seepage of gas and water became more tedious, resulting in decreased relative permeability. Compared with different hydrate pore habits, pore-filling thins the pores, restricting the gas flow than the grain-coating. However, hydrate pore habit barely affects water relative permeability.


---

## Body

## 1. Introduction

Natural gas hydrate (NGH) is an ice-like, crystalline solid in which gas molecules are trapped within water molecules at low temperatures and high pressures [1, 2]. Naturally existing hydrate-bearing sediments (HBS) are widely distributed in continental margins and permafrost [3]. It has been estimated that the total carbon amount stored in the hydrates is at the same order of magnitude as all fossil fuels combined [4]. Therefore, NGH is viewed as a potential energy resource in the future. However, hydrate production is a complex process involving phase transition, multiphase flow, mass transfer, heat transfer, and reservoir deformation [5]. Thus, investigations into safe and efficient hydrate production from HBS have been widely carried out.In-site dissociation is proposed as a feasible and effective method for hydrate production. First, NGH is dissociated by depressurization, thermal stimulation, and inhibitor injection, after which hydrate-dissociated fluids are transported to the surface. During these two processes, HBS permeability is a crucial factor that controls fluid (water or gas) mobility and the ability to transmit pressure and temperature in the reservoir, which affects hydrate productivity [6, 7]. However, HBS permeability varies in hydrate production. Hydrate dissociation releases pore space and changes pore shape, thereby affecting reservoir permeability [8]. Therefore, it is vital to study the permeability change of the reservoir during hydrate exploitation [9]. In addition, in a multiphase flow system, rock fluids interfere with each other due to their different physical and chemical properties at the interface. This interference can be described by relative permeability, related to the saturation of each phase. In other words, relative permeability can better reflect the flowability of fluids in HBS.The research on gas hydrate permeability focuses mainly on the absolute permeability of sediments and the relative permeability of gas and water. The former employs single-phase fluid for assessing absolute permeability varies with hydrate saturation, a condition otherwise called normalized permeability. However, the study of gas/water relative permeability, a multiphase flow system, is more complicated. The existing empirical relative permeability calculation models, such as the Modified Stone model and van Genuchten model, provide only the relationship between the relative permeability of each phase and phase saturation, without considering the secondary influence of NGH dissociation [10, 11]. In addition, the previous experimental studies related to permeability mainly focused on absolute permeability. In contrast, no experiments on the relative permeability of gas and water have been done to our best knowledge [12–15].As an alternative to laboratory experiments, pore network simulation is currently and widely used to study gas/water relative permeability by fitting it with empirical data [16–20]. However, pore network simulation is still insufficient to describe the flow mechanism of multiphase fluids, such as viscous coupling. This model’s simplified porous media structure may also lead to microscopic porosity information loss [21]. Recently, Singh et al. proposed two models that consider and neglect capillary pressure, respectively [22, 23]. These models aid in studying the gas-water multiphase seepage mechanism, but the simplified pore-scale morphology in their derivation may deviate from the actual situation.To better simulate the fluid flow properties and complex pore geometry in HBS, a more accurate pore-scale microsimulation is required. As a mesoscopic simulation method, the Lattice Boltzmann Method (LBM) has been widely used in the past two decades for porous media flow, multicomponent and multiphase flow, chemical reaction, and microscale gas flow [24]. The LBM has a strict microtheoretical basis and some advantages, such as model simplification, ease of handling complex geometry boundary, and good parallelism [25, 26]. However, like laboratory experiments, previous LBM studies of hydrate production have focused on the flow simulation of single-phase fluid; studying the evolution of absolute permeability varies with hydrate saturation [21, 27–30]. Some studies on multiphase flow had considered the effect of capillary number, viscosity ratio, wettability, and fluid distribution on relative permeability [31–33]. So far, the LBM has not been used to study the multiphase flow of gas and water in HBS.The current study used the LBM of multicomponent multiphase (MCMP) pseudopotential model to simulate fluid dual-phase flow in HBS [34]. To effectively simulate fluids’ properties, the Laplace test and contact angle test were used to determine the fluids’ actual density ratio, viscosity ratio, and wettability. Then, a two-dimensional, two-phase Poiseuille flow simulation was applied to verify the model. Finally, the MCMP model was used to study the gas/water dual-phase flow in the reservoir with a particle-filled hydrate for 0.2 saturation. In addition, the evolution characteristics of gas and water relative permeability under different hydrate saturation and hydrate pore habit were discussed to derive the influence of NGH on the multiphase flow.


## 2. Model and Method

### 2.1. Multiphase Flow Simulation

The adopted MCMP Shan-Chen model studies the water and methane immiscible two-phase flow in 2D porous media [35]. The LBM model simulates fluid flow, using the collision and migration of the particle distribution function, fiσx,t for a set of discrete velocities, ei at the lattice site, x, and time, t. The indices, σ, and subscript, i, label the fluid component and discrete velocity direction, respectively. In this study, the Bhatnagar-Gross-Krook model was used, and the evolution equation for σ component is described by:
(1)fiσx+ei∆t,t+∆t−fiσx,t=−1τσfiσx,t−fiσ,eqx,t+Fσ,where τσis the nondimensional relaxation time, related to the kinematic viscosity vσ=cs2τσ−0.5∆t,∆t is the time step (i.e., 1), and Fσ is the forcing term for the σ component. The discrete velocity, ei, along the i direction for D2Q9 lattice model here is written as [36]
(2)e0,e1,e2,e3,e4,e5,e6,e7,e8=c010−101−1−110010−111−1−1.The equilibrium distribution functionfiσ,eqx,t discretized from Maxwell–Boltzmann equilibrium distribution is given by:
(3)fiσ,eq=ωiρσ1+ei∙uσcs2+ei∙uσ22cs4−uσ22cs2,where for the D2Q9 lattice model, weight, ωi, is determined as ω0=4/9,ω1−4=1/9, and ω5−8=1/36. The variable cs is the speed of sound (cs=∆x/∆t/3), while ρσ is the fluid density, obtained as ρσ=∑ifiσ.In the pseudopotential model, the phase separation of the fluid system is due to molecular forces. For the MCMP model, intramolecule force and intermolecule force exist, with the former expressed as:(4)Fσ,σ=−Gσ,σφσxcs2∑i=1Nwei2φσx+eiei,where Gσ,σ is the intramolecule interaction strength and φσx is the effective mass on x site. In our MCMP model, two layers of neighboring nodes are used to calculate the two forces. In one layer of adjacent nodes (Figure 1(a)), only the nearest nodes (nodes 1, 2, 3, and 4) and next-nearest nodes (nodes 5, 6, 7, and 8) are used to calculate the molecular force. This model can lead to the fourth isotropy order (E4). However, two layers of neighboring nodes (Figure 1(b)), with 24 nodes, can reach the eighth isotropy order (E8). Compared with the former, the two layers of neighboring nodes can reduce the spurious currents at the interface and improve model stability [37]. The weight coefficients of the two layers of neighboring nodes are listed in Table 1.Figure 1
Node distributions for calculating the fluid-fluid cohesion force in (a) one-layer and (b) two-layer models [38].
(a)(b)Table 1
Thewei2 of the E8 model used in the current modeling study.
Nodesei2wei21-414/215-824/459-1241/6013-2052/31521-2481/5040The intermolecular force between different fluids can enhance the mutual diffusivity [39]. This force is expressed as:
(5)Fσ,σ¯=−Gσ,σ¯cs2φσx∑i=1Nwei2φσ¯x+eiei,where Gσ,σ¯ is the intermolecule force strength, while Gσσ¯,=Gσ¯,σ=0.5 [40]. Here, we studied two components: σ=1 and σ=2represent the water and methane phases, respectively. Methane is regarded as an ideal fluid; thus, G2,2=0. On the other hand, water is a nonideal fluid, and the Carnahan-Starling (C-S) equation of state is adopted [41]:
(6)P=ρRT1+bρ/4+bρ/42−bρ/431−bρ/43−aρ2,where the critical parameters of the C-S equation of state are as follows:a=0.4963R2Tc2/Pc, and b=0.18727RTc/Pc. We seta=1, b=4, and R=1. To solve the equation of state, we obtained the critical density ρc and critical temperature Tc as 0.13045 and 0.0943, respectively [40]. The effective mass of water component φ1x is given as:
(7)φ1x=2Px−cs2ρxG1,1.When Equation (6) is substituted into Equation (7), the G1,1 cancels out. Hence, the G1,1 value only ensures that square-rooted components of Equation (7) are positive. In addition, the effective mass of methane φ2x is given by:
(8)φ2x=ρ01−expρ2xρ0,where ρ0is a constant coefficient, set as 1, and ρ2x is methane density.To simulate the wettability of each component toward the solid phase, the force of fluids-to-solidFads is calculated. In this study, water and methane were set as wetting and nonwetting fluids, respectively. However, their wettability to the hydrate lack reference data. Thus, we assumed that the fluids have the same wettability to hydrate and solid particles, as follows:
(9)Fads=−Gσ,Scs2φσx∑i=1Nwei2φsx+eiei,whereGσ,S is the controlling parameters for the interface strength between each fluid and the solid. Also, we set G1,S=−1 and G2,S=1. φsx+ei is the effective mass of solid used to adjust contact angle.The body forceG, a steady external force corresponding to the pressure gradient, is used to drive the fluid flow. The effect of gravity on fluid flow was not considered here. Thus, all forces acting on each component can be given as:
(10)Fσ=Fσ,σ+Fσ,k+Fads+G.Here, the exact difference method (EDM) leads to relaxation time independence and has a relatively wide temperature range [42]. The EDM scheme can implement all forces into the LBM framework for a force term in Equation (1). The force term in the evolution equation is expressed as [43]:
(11)Fiσ=fiσ,eqρσ,ueq+∆uσ−fiσ,eqρσ,ueq,where ueq and ∆uσ and can be written as:
(12)ueq=∑σρσuσ∑σρσ,∆uσ=Fσ∆tρσ.In the simulation, the velocity of each component and pressure can be calculated by Equations (13) and (14), respectively:
(13)uσ=∑i=1Nfiσeiρσ+Fσ∆t2ρσ,(14)P=cs2∑σρσ+12cs2Gσ,σφσ2+Gσ,σ¯φσφσ¯+Gσ¯,σφσ¯φσ+Gσ¯,σ¯φσ¯2.


### 2.2. Model Parameters

To fit the actual properties of the fluids, their physical parameters (such as density and viscosity) were predetermined. However, these parameters, especially for the gas, are affected by pressure and temperature. In HBS, the pressure is much higher than the standard atmospheric pressure, resulting in higher density and viscosity. For example, the temperature and pressure in the South China Sea are approximately 14°C and 14 MPa, respectively. Assuming the pressure reduction method is used for hydrate exploitation, and the pressure reduces to 5 MPa, the resulting gas and water physical parameters are given in Table2 (the data refer to AP1700 Material Property Calculation Query Platform).Table 2
Methane and water parameters at 14°C and 5 MPa.
Physical parametersWaterMethaneRatioDensityρkg/m31001.537.226.9Dynamic viscosityμPa∙s1.16×10−31.20×10−596.7Kinetic viscosityvm/s21.16×10−63.23×10−73.6Surface tensionσN/m0.0736


### 2.3. Calculation of Relative Permeability

The relative permeability of waterkrw (or gas, krg) is the water (or gas) effective permeability at a given water saturation Sw, normalized by the water (or gas) permeability at 100% water (or gas) saturation. However, in the Modified Stone equation and van Genuchten model, the gas conductivity at the residual water saturation is used for normalization.In this study, fluid (water or gas) permeability at 100% saturation is used as the reference permeability. Equation (15) calculates the relative permeability:
(15)kr,σ=ke,σka,where kr,σ and ke,σ are the relative permeability and effective permeability, respectively, and ka is the absolute permeability for corresponding HBS.Based on Darcy’s law, the respective equations for calculating the absolute permeabilityka and effective permeability ke,σ by the LBM are determined by:
(16)ka=∑i=1Nus,σiμσNG,ke,σ=∑i=1Num,σiμσNG,where G is the body force, μσ is dynamic viscosity for σ component, N is the total number of grid points at the exit, and us,σi and um,σi are the velocity at each grid point at the exit for σ component in the single-phase and multiphase flows, respectively.The permeability conversion equation between lattice units and physical units is given by:(17)kp=klLpLl2,where kp and kl are the permeability in the physical unit and lattice unit, respectively, while Lp and Ll are the length of the model in the physical unit and lattice unit, respectively.Figure2 shows 2D porous media, including hydrate. The hydrate pore habit of pore-filling and hydrate saturation of Sh=0.2 is an example. In this model, the computational domain was a 360×360 lattice system, and the porosity (except hydrate) was0.65. Solids are composed of similar round particles with a radius of 28 lattices.Figure 2
2D distribution of hydrate and solid particles. The hydrate pore habit is pore-filling. The hydrate and solid particles are shown in black and gray, respectively.The periodic boundary conditions were applied in all directions. The gas and water phases were initially distributed randomly in the pores, and the phase separation occurs automatically. After the phase separation reaches a steady state, the driving force is used to drive two-phase fluid flow from left to right, and then, the outflow flux of each fluid is calculated. Finally, the relative permeability of each phase can be determined from the above equations.In the simulation, it is tedious to find a fixed value for the driving forceG because the viscosity of gas and water is widely different. For example, G is larger for gas flow, disobeying Darcy’s law, while the same G value is smaller for water flow and cannot overcome capillary pressure. Thus, we assume that the effect of capillary number on relative permeability can be negligible, and the G value is changeable. In our simulation, the flow flux may change with time because in gas/water flow system, each phase is usually discontinuous and unsteady. Therefore, every simulation required more than 100000 steps to lower the error in flow flux count.


### 2.4. Model Validation

After the model is determined, the fluid-fluid interaction and the fluid wettability onto the solid should be validated by the Laplace test and contact angle test, respectively. Then, the 2D two-phase Poiseuille flow is simulated to validate that the relative permeability meets the analytical solution.


#### 2.4.1. Laplace Test

The Laplace test was conducted by placing a droplet with different initial radiiR at the center of the lattice domain composed of a 150×150 lattice system. Besides, the periodic boundary condition was applied to the surrounding boundary. Initially, the temperature was set as 0.75Tc, and the density of water (gas) was set as 0.34 (0.0001) and 0.0001 (0.012) for the inside and outside the droplet, respectively, equaling the physical density ratio. To correspond with the kinematic viscosity ratio, the liquid phase and gas phase relaxation times were set as τv=1.2 and τg=0.7, respectively. Provided the density and kinetic viscosity ratios are achieved, the dynamic viscosity is also achievable. Figure 3 shows the density distribution of methane and water along the centerline of the droplet when a steady droplet of water resides within the gas. With the pressure difference between the inside and outside of the droplet and the radiiR, the surface tension can be determined by Laplace’s law:
(18)∆P=σR.Figure 3
The density distributions of methane and water along the centerline of a water droplet.This relationship implies that∆P varies linearly to 1/R. The proportionality coefficient σ is the interfacial tension. Figure 4 depicts ∆P as a function of 1/R where the linear fit shows that the surface tension and σ in the model was 0.0174.Figure 4
The relationship of pressure difference∆Pbetween the water and gas phases and the droplet radii R.


#### 2.4.2. Contact Angle Test

The fluids’ spread on the solid surface can inform on their wettability onto the solid. It depends on the static contact angleθ. For θ<90°, the fluid tends to wet the solid surface, i.e., a wetting fluid. On the other hand, when θ>90°, the fluid alienates the solid surfaces, i.e., a nonwetting fluid. The contact angle between the fluids and the surface can be calculated through Young’s equation, which gives the relationship between the interfacial tension γ12 for fluid-fluid and the interfacial tension γs1 and γs2 for fluids-solid. The Young’s equation is given by:
(19)cosθ=γs2−γs1γ12.Recollect that the gas and the water were set as nonwetting fluid and wetting fluid, respectively. In the MCMP Shan-Chen model, the contact angle is controlled by both fluids near the solid-fluid interface. However, the coeffects are weak for the high-density ratio MCMP model, and the high-density fluid mainly controls the contact angle [34]. Like single component multiphase Shan-Chen model, different static contact angles θ can be achieved by changing the effective solid mass φs. To simulate different contact angles, a 140×100 lattice system with a solid wall at the bottom of the domain was considered. Figure 5 shows the contact angle tests for water. Finally, θ=15°was chosen for the final simulation.Figure 5
Simulation of different contact angles by adjusting the effective solid massφs. The gas is shown as blue, the water red, and the solid and interface between gas and water yellow.


#### 2.4.3. Two-Dimensional, Two-Phase Poiseuille Flow

For immiscible, two-phase flows in a 2D channel, the analytical solutions of speed and relative permeability for each component can be solved by Navier–Stokes equations. Typically, the wetting fluid covers and moves along the solid surface, while the nonwetting fluid flows in the channel center (Figure6). The saturation of each fluid can be defined as the occupied width normalized by the entire width of the channel. The relative permeability of the nonwetting fluid in the center of the channel is affected by the viscosity ratio, M (M=μnw/μw, where μnw and μw are the viscosity of the nonwetting fluid and wetting fluid, respectively).Figure 6
Schematic of 2D immiscible two-phase Poiseuille flow.Huang et al. proposed analytical solutions for each fluid relative permeability as follows [44]:
(20)kr,w=12Sw23−Sw,krn,w=1−Sw32M+1−Sw21−32M.To simulate the Poiseuille flow and validate the simulation results, a240×180 lattice network was established. The no-slip boundary condition was used at the top and bottom of the domain, while the periodic boundary condition was used for the left and right sides. A two-phase cocurrent flow was driven by the body force (G=5×10−9). Figure 7 illustrates the results of numerical simulations and analytical solutions, with a good agreement.Figure 7
Comparison between LBM simulation results and analytical solution for gas and water relative permeability,M=0.010.


## 2.1. Multiphase Flow Simulation

The adopted MCMP Shan-Chen model studies the water and methane immiscible two-phase flow in 2D porous media [35]. The LBM model simulates fluid flow, using the collision and migration of the particle distribution function, fiσx,t for a set of discrete velocities, ei at the lattice site, x, and time, t. The indices, σ, and subscript, i, label the fluid component and discrete velocity direction, respectively. In this study, the Bhatnagar-Gross-Krook model was used, and the evolution equation for σ component is described by:
(1)fiσx+ei∆t,t+∆t−fiσx,t=−1τσfiσx,t−fiσ,eqx,t+Fσ,where τσis the nondimensional relaxation time, related to the kinematic viscosity vσ=cs2τσ−0.5∆t,∆t is the time step (i.e., 1), and Fσ is the forcing term for the σ component. The discrete velocity, ei, along the i direction for D2Q9 lattice model here is written as [36]
(2)e0,e1,e2,e3,e4,e5,e6,e7,e8=c010−101−1−110010−111−1−1.The equilibrium distribution functionfiσ,eqx,t discretized from Maxwell–Boltzmann equilibrium distribution is given by:
(3)fiσ,eq=ωiρσ1+ei∙uσcs2+ei∙uσ22cs4−uσ22cs2,where for the D2Q9 lattice model, weight, ωi, is determined as ω0=4/9,ω1−4=1/9, and ω5−8=1/36. The variable cs is the speed of sound (cs=∆x/∆t/3), while ρσ is the fluid density, obtained as ρσ=∑ifiσ.In the pseudopotential model, the phase separation of the fluid system is due to molecular forces. For the MCMP model, intramolecule force and intermolecule force exist, with the former expressed as:(4)Fσ,σ=−Gσ,σφσxcs2∑i=1Nwei2φσx+eiei,where Gσ,σ is the intramolecule interaction strength and φσx is the effective mass on x site. In our MCMP model, two layers of neighboring nodes are used to calculate the two forces. In one layer of adjacent nodes (Figure 1(a)), only the nearest nodes (nodes 1, 2, 3, and 4) and next-nearest nodes (nodes 5, 6, 7, and 8) are used to calculate the molecular force. This model can lead to the fourth isotropy order (E4). However, two layers of neighboring nodes (Figure 1(b)), with 24 nodes, can reach the eighth isotropy order (E8). Compared with the former, the two layers of neighboring nodes can reduce the spurious currents at the interface and improve model stability [37]. The weight coefficients of the two layers of neighboring nodes are listed in Table 1.Figure 1
Node distributions for calculating the fluid-fluid cohesion force in (a) one-layer and (b) two-layer models [38].
(a)(b)Table 1
Thewei2 of the E8 model used in the current modeling study.
Nodesei2wei21-414/215-824/459-1241/6013-2052/31521-2481/5040The intermolecular force between different fluids can enhance the mutual diffusivity [39]. This force is expressed as:
(5)Fσ,σ¯=−Gσ,σ¯cs2φσx∑i=1Nwei2φσ¯x+eiei,where Gσ,σ¯ is the intermolecule force strength, while Gσσ¯,=Gσ¯,σ=0.5 [40]. Here, we studied two components: σ=1 and σ=2represent the water and methane phases, respectively. Methane is regarded as an ideal fluid; thus, G2,2=0. On the other hand, water is a nonideal fluid, and the Carnahan-Starling (C-S) equation of state is adopted [41]:
(6)P=ρRT1+bρ/4+bρ/42−bρ/431−bρ/43−aρ2,where the critical parameters of the C-S equation of state are as follows:a=0.4963R2Tc2/Pc, and b=0.18727RTc/Pc. We seta=1, b=4, and R=1. To solve the equation of state, we obtained the critical density ρc and critical temperature Tc as 0.13045 and 0.0943, respectively [40]. The effective mass of water component φ1x is given as:
(7)φ1x=2Px−cs2ρxG1,1.When Equation (6) is substituted into Equation (7), the G1,1 cancels out. Hence, the G1,1 value only ensures that square-rooted components of Equation (7) are positive. In addition, the effective mass of methane φ2x is given by:
(8)φ2x=ρ01−expρ2xρ0,where ρ0is a constant coefficient, set as 1, and ρ2x is methane density.To simulate the wettability of each component toward the solid phase, the force of fluids-to-solidFads is calculated. In this study, water and methane were set as wetting and nonwetting fluids, respectively. However, their wettability to the hydrate lack reference data. Thus, we assumed that the fluids have the same wettability to hydrate and solid particles, as follows:
(9)Fads=−Gσ,Scs2φσx∑i=1Nwei2φsx+eiei,whereGσ,S is the controlling parameters for the interface strength between each fluid and the solid. Also, we set G1,S=−1 and G2,S=1. φsx+ei is the effective mass of solid used to adjust contact angle.The body forceG, a steady external force corresponding to the pressure gradient, is used to drive the fluid flow. The effect of gravity on fluid flow was not considered here. Thus, all forces acting on each component can be given as:
(10)Fσ=Fσ,σ+Fσ,k+Fads+G.Here, the exact difference method (EDM) leads to relaxation time independence and has a relatively wide temperature range [42]. The EDM scheme can implement all forces into the LBM framework for a force term in Equation (1). The force term in the evolution equation is expressed as [43]:
(11)Fiσ=fiσ,eqρσ,ueq+∆uσ−fiσ,eqρσ,ueq,where ueq and ∆uσ and can be written as:
(12)ueq=∑σρσuσ∑σρσ,∆uσ=Fσ∆tρσ.In the simulation, the velocity of each component and pressure can be calculated by Equations (13) and (14), respectively:
(13)uσ=∑i=1Nfiσeiρσ+Fσ∆t2ρσ,(14)P=cs2∑σρσ+12cs2Gσ,σφσ2+Gσ,σ¯φσφσ¯+Gσ¯,σφσ¯φσ+Gσ¯,σ¯φσ¯2.


## 2.2. Model Parameters

To fit the actual properties of the fluids, their physical parameters (such as density and viscosity) were predetermined. However, these parameters, especially for the gas, are affected by pressure and temperature. In HBS, the pressure is much higher than the standard atmospheric pressure, resulting in higher density and viscosity. For example, the temperature and pressure in the South China Sea are approximately 14°C and 14 MPa, respectively. Assuming the pressure reduction method is used for hydrate exploitation, and the pressure reduces to 5 MPa, the resulting gas and water physical parameters are given in Table2 (the data refer to AP1700 Material Property Calculation Query Platform).Table 2
Methane and water parameters at 14°C and 5 MPa.
Physical parametersWaterMethaneRatioDensityρkg/m31001.537.226.9Dynamic viscosityμPa∙s1.16×10−31.20×10−596.7Kinetic viscosityvm/s21.16×10−63.23×10−73.6Surface tensionσN/m0.0736


## 2.3. Calculation of Relative Permeability

The relative permeability of waterkrw (or gas, krg) is the water (or gas) effective permeability at a given water saturation Sw, normalized by the water (or gas) permeability at 100% water (or gas) saturation. However, in the Modified Stone equation and van Genuchten model, the gas conductivity at the residual water saturation is used for normalization.In this study, fluid (water or gas) permeability at 100% saturation is used as the reference permeability. Equation (15) calculates the relative permeability:
(15)kr,σ=ke,σka,where kr,σ and ke,σ are the relative permeability and effective permeability, respectively, and ka is the absolute permeability for corresponding HBS.Based on Darcy’s law, the respective equations for calculating the absolute permeabilityka and effective permeability ke,σ by the LBM are determined by:
(16)ka=∑i=1Nus,σiμσNG,ke,σ=∑i=1Num,σiμσNG,where G is the body force, μσ is dynamic viscosity for σ component, N is the total number of grid points at the exit, and us,σi and um,σi are the velocity at each grid point at the exit for σ component in the single-phase and multiphase flows, respectively.The permeability conversion equation between lattice units and physical units is given by:(17)kp=klLpLl2,where kp and kl are the permeability in the physical unit and lattice unit, respectively, while Lp and Ll are the length of the model in the physical unit and lattice unit, respectively.Figure2 shows 2D porous media, including hydrate. The hydrate pore habit of pore-filling and hydrate saturation of Sh=0.2 is an example. In this model, the computational domain was a 360×360 lattice system, and the porosity (except hydrate) was0.65. Solids are composed of similar round particles with a radius of 28 lattices.Figure 2
2D distribution of hydrate and solid particles. The hydrate pore habit is pore-filling. The hydrate and solid particles are shown in black and gray, respectively.The periodic boundary conditions were applied in all directions. The gas and water phases were initially distributed randomly in the pores, and the phase separation occurs automatically. After the phase separation reaches a steady state, the driving force is used to drive two-phase fluid flow from left to right, and then, the outflow flux of each fluid is calculated. Finally, the relative permeability of each phase can be determined from the above equations.In the simulation, it is tedious to find a fixed value for the driving forceG because the viscosity of gas and water is widely different. For example, G is larger for gas flow, disobeying Darcy’s law, while the same G value is smaller for water flow and cannot overcome capillary pressure. Thus, we assume that the effect of capillary number on relative permeability can be negligible, and the G value is changeable. In our simulation, the flow flux may change with time because in gas/water flow system, each phase is usually discontinuous and unsteady. Therefore, every simulation required more than 100000 steps to lower the error in flow flux count.


## 2.4. Model Validation

After the model is determined, the fluid-fluid interaction and the fluid wettability onto the solid should be validated by the Laplace test and contact angle test, respectively. Then, the 2D two-phase Poiseuille flow is simulated to validate that the relative permeability meets the analytical solution.


### 2.4.1. Laplace Test

The Laplace test was conducted by placing a droplet with different initial radiiR at the center of the lattice domain composed of a 150×150 lattice system. Besides, the periodic boundary condition was applied to the surrounding boundary. Initially, the temperature was set as 0.75Tc, and the density of water (gas) was set as 0.34 (0.0001) and 0.0001 (0.012) for the inside and outside the droplet, respectively, equaling the physical density ratio. To correspond with the kinematic viscosity ratio, the liquid phase and gas phase relaxation times were set as τv=1.2 and τg=0.7, respectively. Provided the density and kinetic viscosity ratios are achieved, the dynamic viscosity is also achievable. Figure 3 shows the density distribution of methane and water along the centerline of the droplet when a steady droplet of water resides within the gas. With the pressure difference between the inside and outside of the droplet and the radiiR, the surface tension can be determined by Laplace’s law:
(18)∆P=σR.Figure 3
The density distributions of methane and water along the centerline of a water droplet.This relationship implies that∆P varies linearly to 1/R. The proportionality coefficient σ is the interfacial tension. Figure 4 depicts ∆P as a function of 1/R where the linear fit shows that the surface tension and σ in the model was 0.0174.Figure 4
The relationship of pressure difference∆Pbetween the water and gas phases and the droplet radii R.


### 2.4.2. Contact Angle Test

The fluids’ spread on the solid surface can inform on their wettability onto the solid. It depends on the static contact angleθ. For θ<90°, the fluid tends to wet the solid surface, i.e., a wetting fluid. On the other hand, when θ>90°, the fluid alienates the solid surfaces, i.e., a nonwetting fluid. The contact angle between the fluids and the surface can be calculated through Young’s equation, which gives the relationship between the interfacial tension γ12 for fluid-fluid and the interfacial tension γs1 and γs2 for fluids-solid. The Young’s equation is given by:
(19)cosθ=γs2−γs1γ12.Recollect that the gas and the water were set as nonwetting fluid and wetting fluid, respectively. In the MCMP Shan-Chen model, the contact angle is controlled by both fluids near the solid-fluid interface. However, the coeffects are weak for the high-density ratio MCMP model, and the high-density fluid mainly controls the contact angle [34]. Like single component multiphase Shan-Chen model, different static contact angles θ can be achieved by changing the effective solid mass φs. To simulate different contact angles, a 140×100 lattice system with a solid wall at the bottom of the domain was considered. Figure 5 shows the contact angle tests for water. Finally, θ=15°was chosen for the final simulation.Figure 5
Simulation of different contact angles by adjusting the effective solid massφs. The gas is shown as blue, the water red, and the solid and interface between gas and water yellow.


### 2.4.3. Two-Dimensional, Two-Phase Poiseuille Flow

For immiscible, two-phase flows in a 2D channel, the analytical solutions of speed and relative permeability for each component can be solved by Navier–Stokes equations. Typically, the wetting fluid covers and moves along the solid surface, while the nonwetting fluid flows in the channel center (Figure6). The saturation of each fluid can be defined as the occupied width normalized by the entire width of the channel. The relative permeability of the nonwetting fluid in the center of the channel is affected by the viscosity ratio, M (M=μnw/μw, where μnw and μw are the viscosity of the nonwetting fluid and wetting fluid, respectively).Figure 6
Schematic of 2D immiscible two-phase Poiseuille flow.Huang et al. proposed analytical solutions for each fluid relative permeability as follows [44]:
(20)kr,w=12Sw23−Sw,krn,w=1−Sw32M+1−Sw21−32M.To simulate the Poiseuille flow and validate the simulation results, a240×180 lattice network was established. The no-slip boundary condition was used at the top and bottom of the domain, while the periodic boundary condition was used for the left and right sides. A two-phase cocurrent flow was driven by the body force (G=5×10−9). Figure 7 illustrates the results of numerical simulations and analytical solutions, with a good agreement.Figure 7
Comparison between LBM simulation results and analytical solution for gas and water relative permeability,M=0.010.


## 2.4.1. Laplace Test

The Laplace test was conducted by placing a droplet with different initial radiiR at the center of the lattice domain composed of a 150×150 lattice system. Besides, the periodic boundary condition was applied to the surrounding boundary. Initially, the temperature was set as 0.75Tc, and the density of water (gas) was set as 0.34 (0.0001) and 0.0001 (0.012) for the inside and outside the droplet, respectively, equaling the physical density ratio. To correspond with the kinematic viscosity ratio, the liquid phase and gas phase relaxation times were set as τv=1.2 and τg=0.7, respectively. Provided the density and kinetic viscosity ratios are achieved, the dynamic viscosity is also achievable. Figure 3 shows the density distribution of methane and water along the centerline of the droplet when a steady droplet of water resides within the gas. With the pressure difference between the inside and outside of the droplet and the radiiR, the surface tension can be determined by Laplace’s law:
(18)∆P=σR.Figure 3
The density distributions of methane and water along the centerline of a water droplet.This relationship implies that∆P varies linearly to 1/R. The proportionality coefficient σ is the interfacial tension. Figure 4 depicts ∆P as a function of 1/R where the linear fit shows that the surface tension and σ in the model was 0.0174.Figure 4
The relationship of pressure difference∆Pbetween the water and gas phases and the droplet radii R.


## 2.4.2. Contact Angle Test

The fluids’ spread on the solid surface can inform on their wettability onto the solid. It depends on the static contact angleθ. For θ<90°, the fluid tends to wet the solid surface, i.e., a wetting fluid. On the other hand, when θ>90°, the fluid alienates the solid surfaces, i.e., a nonwetting fluid. The contact angle between the fluids and the surface can be calculated through Young’s equation, which gives the relationship between the interfacial tension γ12 for fluid-fluid and the interfacial tension γs1 and γs2 for fluids-solid. The Young’s equation is given by:
(19)cosθ=γs2−γs1γ12.Recollect that the gas and the water were set as nonwetting fluid and wetting fluid, respectively. In the MCMP Shan-Chen model, the contact angle is controlled by both fluids near the solid-fluid interface. However, the coeffects are weak for the high-density ratio MCMP model, and the high-density fluid mainly controls the contact angle [34]. Like single component multiphase Shan-Chen model, different static contact angles θ can be achieved by changing the effective solid mass φs. To simulate different contact angles, a 140×100 lattice system with a solid wall at the bottom of the domain was considered. Figure 5 shows the contact angle tests for water. Finally, θ=15°was chosen for the final simulation.Figure 5
Simulation of different contact angles by adjusting the effective solid massφs. The gas is shown as blue, the water red, and the solid and interface between gas and water yellow.


## 2.4.3. Two-Dimensional, Two-Phase Poiseuille Flow

For immiscible, two-phase flows in a 2D channel, the analytical solutions of speed and relative permeability for each component can be solved by Navier–Stokes equations. Typically, the wetting fluid covers and moves along the solid surface, while the nonwetting fluid flows in the channel center (Figure6). The saturation of each fluid can be defined as the occupied width normalized by the entire width of the channel. The relative permeability of the nonwetting fluid in the center of the channel is affected by the viscosity ratio, M (M=μnw/μw, where μnw and μw are the viscosity of the nonwetting fluid and wetting fluid, respectively).Figure 6
Schematic of 2D immiscible two-phase Poiseuille flow.Huang et al. proposed analytical solutions for each fluid relative permeability as follows [44]:
(20)kr,w=12Sw23−Sw,krn,w=1−Sw32M+1−Sw21−32M.To simulate the Poiseuille flow and validate the simulation results, a240×180 lattice network was established. The no-slip boundary condition was used at the top and bottom of the domain, while the periodic boundary condition was used for the left and right sides. A two-phase cocurrent flow was driven by the body force (G=5×10−9). Figure 7 illustrates the results of numerical simulations and analytical solutions, with a good agreement.Figure 7
Comparison between LBM simulation results and analytical solution for gas and water relative permeability,M=0.010.


## 3. Simulation Results

### 3.1. Gas and Water Phase Distribution Patterns in HBS

Figure8 depicts gas-water phase distribution patterns in HBS as the water saturationSw changes. Under low water saturation conditions, the water phase is discontinuous, covering the surface of the solid particles and hydrate due to the wettability (Figures 8(a) and 8(b)). However, the gas phase is continuous and flows through the sediments. As the water phase saturation increases (Figure 8(c)), a minute volume of water gradually flows through the sediments. Moreover, the continuous gas phase tends to break up. In our simulation, the water saturation of 0.45 is approximately the dividing point between the continuous and discontinuous gas phases. At Sw=0.56, 0.68, and0.84 (Figures 4(d)–4(f)), the water gradually flows continuously while the gas is disconnected, forming bubbles. Compared with the single gas-phase flow with the same driving force, the gas bubbles flow slowly because of the blockage caused by the relatively higher viscosity water.Figure 8
Gas-water dual-phase distribution patterns in HBS. The gas is shown as blue, water as light green, solid particles as orange, and hydrate as red.Sh=0.2, pore-filling.
(a)
Sw=0.12(b)
Sw=0.3(c)
Sw=0.42(d)
Sw=0.56(e)
Sw=0.68(f)
Sw=0.84


### 3.2. Gas and Water Relative Permeability

Figure9 shows the relative permeability of gas krg and water krw as a function of water saturation Sw. The red and green lines represent krg and krw, respectively. When the water phase saturation is <0.2, the water relative permeability krw is ≈0. At this time, water is bound to the solid surface and cannot flow. Therefore, it can be considered thatSw=0.2 is the residual water saturation under the conditions of our model. When 0.2<Sw<0.45, water starts to flow, although slowly due to lack of a continuous flow path. However, when Sw>0.45, the flow velocity increases significantly as a continuous flow path gradually forms. Finally, the krw approaches 1 until water fills the entire pore space.Figure 9
Relative permeability for gaskrg and water krw as a function of water saturation Sw.Unlike water relative permeability, thekrg decreases with an increase in water saturation. Moreover, when the water saturation is <0.45, the decay rate of the gas relative permeability becomes higher than when the water saturation is >0.45. The following points explain this phenomenon: (a) When the water phase saturation is below the dividing point, the residual water is absorbed on the solid surface and remains immobile. The remaining part of the water flows gently due to the lack of a continuous flow path. In addition, water viscosity is much higher than that of a gas. Due to these reasons, the flow velocity of water is much slower than that of gas. Thus, for the gas phase, increasing water saturation means occupying more pore space. Therefore, the water phase is nearly stagnant and impedes (similar to solid particles) the gas flow. (b) As the water phase saturation exceeds the dividing point, the gas gradually becomes bubbles mixed in the water flow. Here, the gas flow rate gradually matches with the water flow rate. Meanwhile, the water flow has a relatively higher velocity at Sw>0.45, which offsets some negative effects of increasing water saturation on the gas flow rate. Therefore, when the water saturation is >0.45, the gas relative permeability gradually decreases, changing almost linearly with water saturation.


### 3.3. The Main Influence Factors of Gas/Water Relative Permeability

From the gas/water phase distributions and their relative permeability, we observed that high viscosity water consistently and strongly obstructed the gas flow during the entire gas/water multiphase flow. In addition, capillary force is also not conducive to the gas flow. On the contrary, low viscosity gas cannot have a significant inhibitory effect on water flow. The water flow is mainly controlled by the coadsorption of solid particles and hydrates.TheSw at the dividing point is crucial for gas/water relative permeability. For example, the gas phase is continuous at the left side of the dividing point, and the flow velocity is high. However, the water is coadsorbed by the solid particles and hydrate, thereby hindering its movement. Despite that, part of the water beyond the residual water saturation could move; the movement is slow because there is no continuous flow path. On the contrary, gas forms bubbles blocked by water molecules at the right side of the dividing point. However, water’s continuous flow path gradually appears and causes the water flow speed to improve rapidly. This dividing point is for gas; it does not mean water forms a continuous flow path immediately when its saturation is beyond the gas dividing point. The dividing point for water to have a continuous flow path is hard to determine. Therefore, if Sw for the dividing point changes, the relative permeability of gas and water would change.Near the dividing point, the gas relative permeability is always <0.01 because of the 96 times difference in the viscosity of gas and water in our model. As mentioned earlier, the gas changes into bubbles mixed in the water, and its flow rate gradually matches that of the water. However, the water flow rate is 96 times different from that of the gas (single-phase) flow rate under the same driving force. Therefore, the gas relative permeability is always <0.01 if water saturation exceeds the dividing point.


## 3.1. Gas and Water Phase Distribution Patterns in HBS

Figure8 depicts gas-water phase distribution patterns in HBS as the water saturationSw changes. Under low water saturation conditions, the water phase is discontinuous, covering the surface of the solid particles and hydrate due to the wettability (Figures 8(a) and 8(b)). However, the gas phase is continuous and flows through the sediments. As the water phase saturation increases (Figure 8(c)), a minute volume of water gradually flows through the sediments. Moreover, the continuous gas phase tends to break up. In our simulation, the water saturation of 0.45 is approximately the dividing point between the continuous and discontinuous gas phases. At Sw=0.56, 0.68, and0.84 (Figures 4(d)–4(f)), the water gradually flows continuously while the gas is disconnected, forming bubbles. Compared with the single gas-phase flow with the same driving force, the gas bubbles flow slowly because of the blockage caused by the relatively higher viscosity water.Figure 8
Gas-water dual-phase distribution patterns in HBS. The gas is shown as blue, water as light green, solid particles as orange, and hydrate as red.Sh=0.2, pore-filling.
(a)
Sw=0.12(b)
Sw=0.3(c)
Sw=0.42(d)
Sw=0.56(e)
Sw=0.68(f)
Sw=0.84


## 3.2. Gas and Water Relative Permeability

Figure9 shows the relative permeability of gas krg and water krw as a function of water saturation Sw. The red and green lines represent krg and krw, respectively. When the water phase saturation is <0.2, the water relative permeability krw is ≈0. At this time, water is bound to the solid surface and cannot flow. Therefore, it can be considered thatSw=0.2 is the residual water saturation under the conditions of our model. When 0.2<Sw<0.45, water starts to flow, although slowly due to lack of a continuous flow path. However, when Sw>0.45, the flow velocity increases significantly as a continuous flow path gradually forms. Finally, the krw approaches 1 until water fills the entire pore space.Figure 9
Relative permeability for gaskrg and water krw as a function of water saturation Sw.Unlike water relative permeability, thekrg decreases with an increase in water saturation. Moreover, when the water saturation is <0.45, the decay rate of the gas relative permeability becomes higher than when the water saturation is >0.45. The following points explain this phenomenon: (a) When the water phase saturation is below the dividing point, the residual water is absorbed on the solid surface and remains immobile. The remaining part of the water flows gently due to the lack of a continuous flow path. In addition, water viscosity is much higher than that of a gas. Due to these reasons, the flow velocity of water is much slower than that of gas. Thus, for the gas phase, increasing water saturation means occupying more pore space. Therefore, the water phase is nearly stagnant and impedes (similar to solid particles) the gas flow. (b) As the water phase saturation exceeds the dividing point, the gas gradually becomes bubbles mixed in the water flow. Here, the gas flow rate gradually matches with the water flow rate. Meanwhile, the water flow has a relatively higher velocity at Sw>0.45, which offsets some negative effects of increasing water saturation on the gas flow rate. Therefore, when the water saturation is >0.45, the gas relative permeability gradually decreases, changing almost linearly with water saturation.


## 3.3. The Main Influence Factors of Gas/Water Relative Permeability

From the gas/water phase distributions and their relative permeability, we observed that high viscosity water consistently and strongly obstructed the gas flow during the entire gas/water multiphase flow. In addition, capillary force is also not conducive to the gas flow. On the contrary, low viscosity gas cannot have a significant inhibitory effect on water flow. The water flow is mainly controlled by the coadsorption of solid particles and hydrates.TheSw at the dividing point is crucial for gas/water relative permeability. For example, the gas phase is continuous at the left side of the dividing point, and the flow velocity is high. However, the water is coadsorbed by the solid particles and hydrate, thereby hindering its movement. Despite that, part of the water beyond the residual water saturation could move; the movement is slow because there is no continuous flow path. On the contrary, gas forms bubbles blocked by water molecules at the right side of the dividing point. However, water’s continuous flow path gradually appears and causes the water flow speed to improve rapidly. This dividing point is for gas; it does not mean water forms a continuous flow path immediately when its saturation is beyond the gas dividing point. The dividing point for water to have a continuous flow path is hard to determine. Therefore, if Sw for the dividing point changes, the relative permeability of gas and water would change.Near the dividing point, the gas relative permeability is always <0.01 because of the 96 times difference in the viscosity of gas and water in our model. As mentioned earlier, the gas changes into bubbles mixed in the water, and its flow rate gradually matches that of the water. However, the water flow rate is 96 times different from that of the gas (single-phase) flow rate under the same driving force. Therefore, the gas relative permeability is always <0.01 if water saturation exceeds the dividing point.


## 4. Discussion

The simulation results show the evolution of gas/water relative permeability, and their phase distributions vary with water saturation under a certain hydrate saturation and hydrate pore habit. However, the hydrate dissociation may bring changes to the microscopic pore structure of sediments, which affects the fluid microseepage behavior and their relative permeability evolution. To analyze the influence of hydrate growth behavior (i.e., saturation and pore habit) on gas and water relative permeability, we simulated the multiphase flow behaviors atSh=0.2, 0.4, and 0.6 and the hydrate pore habit of grain-coating and pore-filling.


### 4.1. Effect of Hydrate Saturation on Gas/Water Relative Permeability

Figure10 shows gas-water dual-phase distribution patterns in HBS with 0.4 and 0.6 hydrate saturation for the hydrate pore habit of pore-filling, respectively. As the hydrate saturation increases, the fluid distribution patterns change in the following aspects: (a) Gas phase is more likely to exist as bubbles and challenging to overcome capillary resistance. Moreover, the gas phase is more likely to change from continuous to discontinuous as hydrate saturation increases under the same water saturation conditions. (b) The dividing point between the continuous and discontinuous gas phases changes from 0.45 to 0.30 as the hydrate saturation increases from 0.2 to 0.6 (Figure 10(a)). (c) More water molecules are bound to the solid surface, thereby increasing residual water saturation. However, once a continuous flow path has been established for the water phase, the flow velocity increases rapidly and is less affected by external factors. (d) The gas bubbles which could flow gradually become residual gas as hydrate saturation increases (Figure 10(c)). In general, the increased hydrate saturation hinders the fluid flow. However, water flow is less affected by hydrate saturation than gas flow.Figure 10
Gas and water dual-phase distribution patterns atSh=0.4 and 0.6 for the hydrate pore habit of pore-filling. The gas is shown as blue, water as light green, solid particles as orange, and hydrate as red.
(a)
Sh=0.4(b)
Sh=0.6The evolution of gas/water relative permeability caused by different hydrate saturation is shown in Figure11. The black spot is the pore network simulation reported by Mahabadi et al. [17], the green spot is the analysis model studied by Singh et al. [23], and the solid red line is from the current research using the LBM. We observed that the simulation result of the present study matches with the analysis model but with a slight deviation from the pore network model. In pore network simulation, the relative permeability, especially gas relative permeability, is higher than that of the other two methods. It may be because the pore network simulation is difficult to capture the viscous interaction between different fluids, ignoring the obstructive effect of water on gas.Figure 11
Relative permeability for gaskrg and water krw as a function of water saturation Sw for the hydrate pore habit of pore-filling.Furthermore, we observed (Figure12) that the gas relative permeability (krg) reduces as hydrate saturation increases. As the hydrate saturation increases from 0.2 to 0.6, the gas relative permeability trend shrinks toward the left. In these figures, krg=0.01 is an important point that corresponds to water saturationSw near the dividing point (Section 3.3). However, the corresponding water saturationSw for the dividing point reduces as the hydrate saturation increases. Thus, the declining rate of the gas relative permeability curves is significantly faster. Besides, as hydrate saturation increases, the gas bubbles which could flow gradually become residual gas, making the gas relative permeability reaches zero earlier. As Figure 11 shows, when hydrate saturation increases from 0.2 to 0.6, the residual gas increases from 0 to 0.2.Figure 12
Gas and water dual-phase distribution patterns inSh=0.2,0.4, and 0.6 for the hydrate pore habit of grain-coating. Water is shown as light green, solid particles as orange, and hydrate as red.
(a)
Sh=0.2(b)
Sh=0.4(c)
Sh=0.6For the relative permeability of water, as hydrate saturation increases from 0.2 to 0.6, the water relative permeability trend shrank to the right. Water saturation is required for the initial water flow to increase, and the residual water increases from 0.2 to 0.3 when the hydrate saturation increases from 0.2 to 0.6 (Figure11). Besides, the water relative permeability increases slowly initially. In our simulation, the effect of hydrate on water relative permeability is mainly reflected in these low water saturation conditions. When Sw>0.6, the curves of water relative permeability do not change significantly with hydrate saturation. In other words, hydrate does not greatly impact water relative permeability when water saturation is high.


### 4.2. Effect of Hydrate Pore Habit on Gas/Water Relative Permeability

This section analyzes the effect of hydrate pore habit on gas/water relative permeability evolution behaviors.The two-phase distribution patterns for hydrate pore habit of grain-coating are shown in Figure12. Compared with Figure 10, the hydrate pore habit of grain-coating also encounters obstacles on the fluid flow, like pore-filling. However, the extent of influence on gas/water seepage characteristics between these two habits is unclear from these phase distribution patterns.Figure13 shows the comparison of gas and water relative permeability in their hydrate pore habits. The red and blue lines show the pore habits of pore-filling and grain-coating, respectively. Except for Sh=0.4 and Sw>0.5, the gas relative permeability for grain-coating is almost higher than that of pore-filling. Because the hydrate occupied the center of pores in the pore-filling, the pores shrink than that of grain-coating, and a more prominent obstacle to the gas flow occurred. The case of Sh=0.4 and Sw>0.5 appearing is an abnormal situation. This scenario is due to the pore habit of grain-coating, forming larger pores than that of pore-filling, which causes gas bubbles to become residual gas more easily. However, the gas relative permeability for the hydrate pore habit of grain-coating is generally higher than that of pore-filling. However, for water relative permeability, the trends of the hydrate types are similar to hydrate saturation, increasing from 0.2 to 0.6. Therefore, the hydrate pore habit has the same effect on water relative permeability.Figure 13
Comparison of gas/water relative permeability for hydrate pore habit of grain-coating and pore-filling.The hydrate pore habit has the same effect on water relative permeability because the water seepage can only be affected by hydrates when water saturation is low. Regardless, the pore-filling and grain-coating have a strong adsorption effect on water flow in low water saturation conditions. Moreover, the two types of hydrate pore habits have little difference in the residual water saturation and the water saturation required to form a continuous flow path. In high water saturation conditions, as long as the continuous flow path for water has been formed, hydrate pore habit and hydrate saturation cannot significantly influence the water relative permeability. On the contrary, gas relative permeability is mainly influenced by viscous coupling between gas and water, capillary resistance, and residual gas saturation. Although pore-filling cannot induce more residual gas, it, however, causes pore thinning. Therefore, this form of pore-filling causes more immense capillary pressure and makes more gas form discrete bubbles. For these reasons, gas relative permeability for the pore habit of pore-filling is always lower than the pore habit of grain-coating.


## 4.1. Effect of Hydrate Saturation on Gas/Water Relative Permeability

Figure10 shows gas-water dual-phase distribution patterns in HBS with 0.4 and 0.6 hydrate saturation for the hydrate pore habit of pore-filling, respectively. As the hydrate saturation increases, the fluid distribution patterns change in the following aspects: (a) Gas phase is more likely to exist as bubbles and challenging to overcome capillary resistance. Moreover, the gas phase is more likely to change from continuous to discontinuous as hydrate saturation increases under the same water saturation conditions. (b) The dividing point between the continuous and discontinuous gas phases changes from 0.45 to 0.30 as the hydrate saturation increases from 0.2 to 0.6 (Figure 10(a)). (c) More water molecules are bound to the solid surface, thereby increasing residual water saturation. However, once a continuous flow path has been established for the water phase, the flow velocity increases rapidly and is less affected by external factors. (d) The gas bubbles which could flow gradually become residual gas as hydrate saturation increases (Figure 10(c)). In general, the increased hydrate saturation hinders the fluid flow. However, water flow is less affected by hydrate saturation than gas flow.Figure 10
Gas and water dual-phase distribution patterns atSh=0.4 and 0.6 for the hydrate pore habit of pore-filling. The gas is shown as blue, water as light green, solid particles as orange, and hydrate as red.
(a)
Sh=0.4(b)
Sh=0.6The evolution of gas/water relative permeability caused by different hydrate saturation is shown in Figure11. The black spot is the pore network simulation reported by Mahabadi et al. [17], the green spot is the analysis model studied by Singh et al. [23], and the solid red line is from the current research using the LBM. We observed that the simulation result of the present study matches with the analysis model but with a slight deviation from the pore network model. In pore network simulation, the relative permeability, especially gas relative permeability, is higher than that of the other two methods. It may be because the pore network simulation is difficult to capture the viscous interaction between different fluids, ignoring the obstructive effect of water on gas.Figure 11
Relative permeability for gaskrg and water krw as a function of water saturation Sw for the hydrate pore habit of pore-filling.Furthermore, we observed (Figure12) that the gas relative permeability (krg) reduces as hydrate saturation increases. As the hydrate saturation increases from 0.2 to 0.6, the gas relative permeability trend shrinks toward the left. In these figures, krg=0.01 is an important point that corresponds to water saturationSw near the dividing point (Section 3.3). However, the corresponding water saturationSw for the dividing point reduces as the hydrate saturation increases. Thus, the declining rate of the gas relative permeability curves is significantly faster. Besides, as hydrate saturation increases, the gas bubbles which could flow gradually become residual gas, making the gas relative permeability reaches zero earlier. As Figure 11 shows, when hydrate saturation increases from 0.2 to 0.6, the residual gas increases from 0 to 0.2.Figure 12
Gas and water dual-phase distribution patterns inSh=0.2,0.4, and 0.6 for the hydrate pore habit of grain-coating. Water is shown as light green, solid particles as orange, and hydrate as red.
(a)
Sh=0.2(b)
Sh=0.4(c)
Sh=0.6For the relative permeability of water, as hydrate saturation increases from 0.2 to 0.6, the water relative permeability trend shrank to the right. Water saturation is required for the initial water flow to increase, and the residual water increases from 0.2 to 0.3 when the hydrate saturation increases from 0.2 to 0.6 (Figure11). Besides, the water relative permeability increases slowly initially. In our simulation, the effect of hydrate on water relative permeability is mainly reflected in these low water saturation conditions. When Sw>0.6, the curves of water relative permeability do not change significantly with hydrate saturation. In other words, hydrate does not greatly impact water relative permeability when water saturation is high.


## 4.2. Effect of Hydrate Pore Habit on Gas/Water Relative Permeability

This section analyzes the effect of hydrate pore habit on gas/water relative permeability evolution behaviors.The two-phase distribution patterns for hydrate pore habit of grain-coating are shown in Figure12. Compared with Figure 10, the hydrate pore habit of grain-coating also encounters obstacles on the fluid flow, like pore-filling. However, the extent of influence on gas/water seepage characteristics between these two habits is unclear from these phase distribution patterns.Figure13 shows the comparison of gas and water relative permeability in their hydrate pore habits. The red and blue lines show the pore habits of pore-filling and grain-coating, respectively. Except for Sh=0.4 and Sw>0.5, the gas relative permeability for grain-coating is almost higher than that of pore-filling. Because the hydrate occupied the center of pores in the pore-filling, the pores shrink than that of grain-coating, and a more prominent obstacle to the gas flow occurred. The case of Sh=0.4 and Sw>0.5 appearing is an abnormal situation. This scenario is due to the pore habit of grain-coating, forming larger pores than that of pore-filling, which causes gas bubbles to become residual gas more easily. However, the gas relative permeability for the hydrate pore habit of grain-coating is generally higher than that of pore-filling. However, for water relative permeability, the trends of the hydrate types are similar to hydrate saturation, increasing from 0.2 to 0.6. Therefore, the hydrate pore habit has the same effect on water relative permeability.Figure 13
Comparison of gas/water relative permeability for hydrate pore habit of grain-coating and pore-filling.The hydrate pore habit has the same effect on water relative permeability because the water seepage can only be affected by hydrates when water saturation is low. Regardless, the pore-filling and grain-coating have a strong adsorption effect on water flow in low water saturation conditions. Moreover, the two types of hydrate pore habits have little difference in the residual water saturation and the water saturation required to form a continuous flow path. In high water saturation conditions, as long as the continuous flow path for water has been formed, hydrate pore habit and hydrate saturation cannot significantly influence the water relative permeability. On the contrary, gas relative permeability is mainly influenced by viscous coupling between gas and water, capillary resistance, and residual gas saturation. Although pore-filling cannot induce more residual gas, it, however, causes pore thinning. Therefore, this form of pore-filling causes more immense capillary pressure and makes more gas form discrete bubbles. For these reasons, gas relative permeability for the pore habit of pore-filling is always lower than the pore habit of grain-coating.


## 5. Conclusion

We used LBM to study the effect of hydrate on gas/water relative permeability. Gas and water multiphase flows for the hydrate saturation ofSh=0.2,0.4, and 0.6, with the hydrate pore habit of grain-coating and pore-filling, are simulated, and their relative permeability is analyzed. The main conclusions are as follows:
(1)
Our research has the same trend as these two methods compared to the pore network simulation and analysis model. This observation shows that the LBM is valid in simulating gas and water multiphase flow in HBS(2)
During the multiphase flow process of gas and water, high viscosity water has a strong, obstructive effect on gas flow due to the viscous coupling, but the gas phase has little effect on water flow. The water flow is mainly controlled by the coadsorption of solid particles and hydrates. Whether the water forms a continuous flow path is another critical factor. In our model, the corresponding water saturationSw of the dividing point is approximately 0.45 when hydrate saturation is 0.20. On the left side of the dividing point, the gas relative permeability quickly decreases from 1 to 0.01. In contrast, the water relative permeability starts as water saturation exceeds 0.2, but it is usually <0.1. On the right side of the dividing point, the gas relative permeability decreases gradually from 0.01 to 0, while the water relative permeability quickly increases from 0.1 to 1(3)
Gas and water multiphase seepage becomes more difficult as hydrate saturation increases, and their relative permeability decreases. There are three main reasons for the decrease in gas relative permeability: (i) The dividing water saturationSw between continuous and discontinuous for the gas phase decreases. In our simulation, as hydrate saturation increases from 0.2 to 0.6, the dividing water saturationSw decreases from 0.45 to 0.30. (ii) Increased hydrate saturation causes pores to narrow, thus resulting in an enormous capillary resistance. (iii) Mobile gas bubbles gradually become residual gas as hydrate saturation increases. As hydrate saturation increases from 0.2 to 0.6, the residual gas saturation increases from 0 to 0.2. Moreover, the reduction of water relative permeability is mainly due to the coadsorption of solid particles and hydrates. The coadsorption increases the water saturation required for initial water flow and the water saturation required to form a continuous flow path. As hydrate saturation increases from 0.2 to 0.6, the residual water saturation increases from 0.2 to 0.3(4)
For different hydrate pore habits, the gas relative permeability for pore-filling is lower than that of grain-coating because hydrate occupied in the center of pores thins the pore further. However, the water relative permeability is less affected by the hydrate pore habit



---

*Source: 1396323-2021-10-25.xml*