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

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 saturation Sw 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.


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][13][14][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][17][18][19][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 gaswater 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][28][29][30]. Some studies on multiphase flow had considered the effect of capillary number, viscosity ratio, wettability, and fluid distribution on relative permeability [31][32][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, twophase 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 particlefilled 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.

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, f σ i ðx, tÞ for a set of discrete velocities, e i 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: where τ σ is the nondimensional relaxation time, related to the kinematic viscosity ðv σ = c 2 s ðτ σ − 0:5ÞΔtÞ,Δt is the time step (i.e., 1), and F σ is the forcing term for the σ component. The discrete velocity, e i , along the i direction for D 2 Q 9 lattice model here is written as [36] e 0 , e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 , e 8 ½ = c The equilibrium distribution function f σ,eq i ðx, tÞ discretized from Maxwell-Boltzmann equilibrium distribution is given by: where for the D 2 Q 9 lattice model, weight, ω i , is determined as ω 0 = 4/9, ω 1−4 = 1/9, and ω 5−8 = 1/36. The variable c s is the speed of sound (c s = ðΔx/ΔtÞ/ ffiffi ffi 3 p ), while ρ σ is the fluid density, obtained as ρ σ = ∑ i f σ i . 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: 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 nextnearest 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.
The intermolecular force between different fluids can enhance the mutual diffusivity [39]. This force is expressed as: where G σ, σ is the intermolecule force strength, while G σ σ, = G σ,σ =0.5 [40]. Here, we studied two components: σ = 1 and σ = 2 represent the water and methane phases, respectively. Methane is regarded as an ideal fluid; thus, G 2,2 = 0. On the other hand, water is a nonideal fluid, and the Carnahan-Starling (C-S) equation of state is adopted [41]: where the critical parameters of the C-S equation of state are as follows: a = 0:4963R 2 T 2 c /P c , and b = 0:18727RT c /P c . We set a = 1, b = 4, and R = 1. To solve the equation of state, we obtained the critical density ρ c and critical temperature T c as 0.13045 and 0.0943, respectively [40]. The effective mass of water component φ 1 ðxÞ is given as: When Equation (6) is substituted into Equation (7), the G 1,1 cancels out. Hence, the G 1,1 value only ensures that square-rooted components of Equation (7) are positive. In addition, the effective mass of methane φ 2 ðxÞ is given by: where ρ 0 is a constant coefficient, set as 1, and ρ 2 ðxÞ is methane density.
To simulate the wettability of each component toward the solid phase, the force of fluids-to-solid F ads 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: where G σ,S is the controlling parameters for the interface strength between each fluid and the solid. Also, we set G 1,S = −1 and G 2,S = 1. φ s ðx + e i Þ is the effective mass of solid used to adjust contact angle. The body force G, 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: 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]: where u eq and Δu σ and can be written as: In the simulation, the velocity of each component and pressure can be calculated by Equations (13) and (14), respectively:

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 Table 2 (the data refer to AP1700 Material Property Calculation Query Platform).

Calculation of Relative
Permeability. The relative permeability of water k rw (or gas, k rg ) is the water (or gas) effective permeability at a given water saturation S w, 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: where k r,σ and k e,σ are the relative permeability and effective permeability, respectively, and k a is the absolute permeability for corresponding HBS. Based on Darcy's law, the respective equations for calculating the absolute permeability ðk a Þ and effective permeability ðk e,σ Þ by the LBM are determined by: where G is the body force, μ σ is dynamic viscosity for σ component, N is the total number of grid points at the exit, and u s,σi and u m,σ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: where k p and k l are the permeability in the physical unit and lattice unit, respectively, while L p and L l are the length of the model in the physical unit and lattice unit, respectively.   Geofluids Figure 2 shows 2D porous media, including hydrate. The hydrate pore habit of pore-filling and hydrate saturation of S h = 0:2 is an example. In this model, the computational domain was a 360 × 360 lattice system, and the porosity (except hydrate) was 0:65. Solids are composed of similar round particles with a radius of 28 lattices.
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 twophase 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 force G 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.

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 radii R 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:75T c , 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 radii R, the surface tension can be determined by Laplace's law: This relationship implies that ΔP varies linearly to 1/R. The proportionality coefficient σ is the interfacial tension.

Geofluids
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.
To simulate the Poiseuille flow and validate the simulation results, a 240 × 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 8 depicts gas-water phase distribution patterns in HBS as the water saturation S w 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 S w = 0:56, 0:68, and 0: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.

Gas and Water
Relative Permeability. Figure 9 shows the relative permeability of gas ðk rg Þ and water ðk rw Þ as a function of water saturation ðS w Þ. The red and green lines represent k rg and k rw , respectively. When the water phase saturation is <0.2, the water relative permeability ðk rw Þ is ≈0. At this time, water is bound to the solid surface and cannot flow. Therefore, it can be considered that S w = 0:2 is the residual water saturation under the conditions of our model. When 0:2 < S w < 0:45, water starts to flow, although slowly due to lack of a continuous flow path. However, when S w > 0:45, the flow velocity increases significantly as a continuous flow path gradually forms. Finally, the k rw approaches 1 until water fills the entire pore space.
Unlike water relative permeability, the k rg 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 S w > 0:45, which offsets some 8 Geofluids 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.

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. The S w 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 S w 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 (singlephase) flow rate under the same driving force. Therefore, the gas relative permeability is always <0.01 if water saturation exceeds the dividing point.

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 at S h = 0:2, 0:4, and 0:6 and the hydrate pore habit of grain-coating and pore-filling.

Effect of Hydrate Saturation on Gas/Water Relative
Permeability. Figure 10 shows gas-water dual-phase distribution patterns in HBS with 0.4 and 0.6 hydrate saturation Pore-filling S h =0.2 Figure 11: Relative permeability for gas ðk rg Þ and water ðk rw Þ as a function of water saturation ðS w Þ for the hydrate pore habit of pore-filling.
9 Geofluids 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.
The evolution of gas/water relative permeability caused by different hydrate saturation is shown in Figure 11. 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.
Furthermore, we observed ( Figure 12) that the gas relative permeability (k rg ) 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 However, the corresponding water saturation S w 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. For 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 ( Figure 11). 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 S w > 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.

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 Figure 12. 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. Figure 13 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 S h = 0:4 and S w > 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 S h = 0:4 and S w > 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.
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.

Conclusion
We used LBM to study the effect of hydrate on gas/water relative permeability. Gas and water multiphase flows for the hydrate saturation of S h = 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 saturation S w 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 saturation S w 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 saturation S w 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

Conflicts of Interest
The authors declare that they have no conflicts of interest.