# Properties and Model of Pore-Scale Methane Displacing Water in Hydrate-Bearing Sediments
## 3 Numerical Model Description  
### Governing Equation  
To visualize the constrained behavior of the two phases within the microchannel, simulation operations have been conducted in COMSOL using the level-set method. Tracking the shape and position of the interface is accomplished by the time evolution of the level-set function. The model incorporates interfacial tension as a mass force into the Navier-Stokes equations, which are utilized to depict mass and momentum transfer in the two-phase methane and water fluids [27]. \[\rho\frac{\partial u}{\partial t}+\rho(u\cdot\
abla)u=\
abla\cdot\left[-pI+ \mu\Big{(}\
abla u+(\
abla u)^{T}\Big{)}\right]+F_{st} \tag{1}\]  
\[\frac{\partial\rho}{\partial t}+\
abla\cdot(\rho u)=0 \tag{2}\]  
where \(\rho\) is expressed as the density of the fluid, kg/cm\({}^{3}\). \(\mu\) represents the dynamic viscosity of the fluid, Pa-s. \(u\) is the speed, m/s. p represents the inlet pressure, Pa. \(I\) is the unit matrix. \(F_{st}\) is the surface tension of the two phases, N/m. The flow velocity at the minimum opening of the simulated process in this paper is small, and the Reynolds number is much less than 1, which can ignore the inertial term of the flow. In this way, the fluid flow problem is transformed into solving the Navier-Stokes equations under specific conditions. The above equations can be rewritten as follows [28]:  
\[\frac{\partial u}{\partial t}+(u\cdot\
abla)u+\frac{1}{\rho}\
abla p-\mu\
abla ^{2}u=F_{st} \tag{3}\]  
\[\rho\
abla\cdot u=0 \tag{4}\]  
where \(p\) is expressed as pressure, Pa. \(u\) is expressed as the flow velocity, m/s. \(\mu\) is denoted as the dynamic viscosity of the fluid, Pa-s. \(\rho\) is expressed as the density of the fluid, kg/m\({}^{3}\)  
Figure 4: Hydrate mining and pore frequency distribution. The respective volume fractions of the gas-water two-phase flow are defined in a separated two-phase flow by the level-set function \(\varphi\), where \(0\leq\varphi\leq 1\). In the two-phase flow simulation process, \(\varphi=0\) indicates the aqueous phase, and \(\varphi=1\) indicates the gas phase. When \(0<\varphi<1\), \(\varphi\) is used to describe a phase interface with a certain thickness. It is a value that presents a gradient change from 0 to 1. A contour surface with \(\varphi=0.5\) is typically considered the phase interface in post-processing. The two-phase flow interface is tracked using the reinitialized level-set function \(\varphi\) values to study the dynamic characteristics of the gas-water two-phase flow interface. The kinematic equation of the level-set method for tracking the interface of the gas-water two-phase flow can be expressed as follows [29]:  
\[\frac{\partial\phi}{\partial\mathrm{t}}+u\cdot\
abla\phi=\gamma\
abla\cdot \left(\varepsilon\
abla\phi-\phi(1-\phi)\frac{\
abla\phi}{|\
abla\phi|}\right) \tag{5}\]  
where \(\gamma\) is the reinitialization parameter; \(\varepsilon\) is the interfacial thickness control parameter; \(\partial\phi/\partial\mathrm{t}\) denotes an accumulation term over time; and \(u\cdot\phi\) denotes a convective term. The right-hand side of the level-set equation represents the diffusion term. The choice of both \(\gamma\) and \(\varepsilon\) affects the astringency of the model. To improve astringency when calculating the level-set equations, \(\varepsilon=h_{c}/2\) is usually used, where \(h_{c}\) is the grid size of the interface region, and \(\gamma\) is the maximum velocity occurring in the model. A uniform mesh size is more conducive to computational convergence, and an appropriate increase in \(\varepsilon\) can accelerate the model. Since the level-set function is a smooth step function, the density (\(\rho\)) and the dynamic viscosity (\(\mu\)) of the two-phase mixing parameters can be determined by the following equation [30]:  
\[ \begin{array}{l}\rho=\rho_{g}+\left(\rho_{w}-\rho_{g}\right)\varphi\\ \mu=\mu_{g}+\left(\mu_{w}-\mu_{g}\right)\varphi\end{array}