Inhibition mechanism of solid phase migration and flow channel jamming on the natural gas hydrate dissociation and gas production

Abstract

In the process of natural gas hydrate exploitation, solid-phase particles tend to migrate in pore channels, resulting in localized flow channel jamming and consequently destabilizing gas production. This presents a significant challenge for the commercial exploitation of hydrates. Therefore, uncovering the influence of solid phase migration in the flow channel holds great importance for achieving more efficient and safe hydrate exploitation. In this study, a gas hydrate exploitation model was established based on the two-fluid model theory, incorporating the formation of ice phase and accounting for solid phase migration-jamming. The results indicate that both gas production rate and gas hydrate dissociation rate decrease due to the influence of hydrate phase and ice phase migration and jamming. This effect is attributed to the increased solid phase saturation as solid particles flow toward the near-wellbore area. Consequently, the absolute permeability of the reservoir is inhibited, limiting the gas–liquid flow capacity, particularly in instances where the flow channel becomes jamming. Convective heat transfer within the reservoir diminishes, leading to a rapid drop in the local reservoir temperature. Simultaneously, the quantity of the ice phase increases within this localized area. However, the migration of solid particles has minimal impact on changes in the reservoir's pressure field. Therefore, variations in the pressure field are not the primary cause of gas production reduction resulting from solid phase migration influence.

Introduction

Natural gas hydrate (NGH) is one of the most promising new energy sources. The process of exploitation of NGH is complex, and it needs to be studied in deep. At present, the generally accepted methods for the exploitation of NGH mainly include: depressurization, heating, inhibitor injection and CO₂ replacement (Choudhary and Phirani, 2023; Shao et al., 2024; Zhang et al., 2024). The field test data and numerical simulation results all proved that the depressurization method is the most effective method of gas hydrate production (Lei et al., 2022; Moridis et al., 2008; Yamamoto et al., 2022; Yin and Linga, 2019). However, the commercial exploitation of hydrates is still not possible, because the production of large amounts of solid particles and water can easily lead to the jamming of mining well and flow channels (Konno et al., 2017); the methane gas leakage and geological hazards need to be predicted and prevented (Boudreau et al., 2015; Yao and Wu, 2021; Zhu et al., 2025); and the economic efficiency cannot be evaluated effectively now. Therefore, further in-depth and comprehensive research is needed on the depressurization process of natural gas hydrates, particularly regarding the mechanism studies of solid particle generation, transport, and the flow channel jamming issues caused by them.

The dissociation of natural gas hydrates involves a significant endothermic reaction, which induces the Joule–Thomson cooling effect. This effect causes a rapid decline in the current temperature state of natural gas hydrates, leading to the formation of ice phase under certain conditions (Mahmood and Guo, 2021, 2023). Ohno et al.'s experimental results showed that the effect of ice on hydrate dissociation rate can be divided into two stages. In the initial stage, the dissociation of hydrate produces gas and ice, first the dissociation rate decreases significantly, and then the dissociation rate enters the rising stage (Abbasi et al., 2022). The inhibitory effect of the ice phase on the dissociation rate of hydrate particles beneath the ice shell has been verified at the molecular scale (Hassanpouryouzband et al., 2020). The impact of the ice phase on hydrate exploitation cannot be disregarded under certain reservoir conditions. In our previous study, a numerical model was proposed to quantitatively describe the effect of ice shell on the dissociation rate of hydrate (Yu et al., 2018b). This model took into account the change in hydrate dissociation surface area due to the covering of ice attached to the surface of hydrate particle.

During hydrates exploitation process, solid particles, such as hydrates, ice, and sand, could migrate with the gas and liquid flow in the flow channel. The productivity of the well is susceptible to disruption due to the emergence of solid-phase materials near the producing wellbore, causing potential flow restrictions or jamming (Xia et al., 2022). Previous studies have found that the permeability of the reservoir decreases with the fluid velocity increasing because the larger fluid velocity increases the drag force on the solid particles and enhances the solid migration (Fang et al., 2017; Gruesbeck and Collins, 1982). Under the constant back pressure, the fluid velocity in the reservoir decreases due to the solid particles migration (Shan and Zhou, 2020). Cai et al. (Cai et al., 2020) verified that the reservoir porosity and permeability both decrease with increasing pressure difference. Moreover, as a special case in the multiphase flow, flow channel jamming could form due to the solid particles migration and deposition (Hauge et al., 2016). When the jamming occurs, the solid particles stop, the gas productive rate and other parameters could be significantly affected. Therefore, the mechanism of jamming cannot be ignored when the solid particles migration is considered in the hydrate exploitation simulation.

The jamming mechanism was not studied and published as a scientific question until the middle of the last century (Wang et al., 2011). To et al. (To et al., 2001) first proposed a gravity-drive jamming rate model by using two-dimensional particles falling into a funnel experiment. Subsequently, through a large number of experiments, it was found that the occurrence of jamming was related to the shape and size of particles, opening size, vibration, and other factors (Magalhães et al., 2010; Mankoc et al., 2009). In recent years, researchers have begun to study the mechanism of fluid-driven jamming (Guariguata et al., 2012; Janda et al., 2012; Lafond, 2014; Lafond et al., 2013; Nan et al., 2018, 2024). Guariguata et al. (Guariguata et al., 2012) analyzed the influence of flow velocity, channel size, particle size, and inlet flow characteristics on the system jamming. According to the two-dimensional empirical formula established by Janda et al. (Janda et al., 2012) for the particle quantity rate of solid particles passing through the throttling hole driven by gravity, Lafond et al. (Lafond, 2014) deduced the three-dimensional calculation formula for the particle quantity rate of particles passing through the throttling hole driven by fluid. These research findings provide a foundation basis for analyzing the effects of solid particle migration and jamming within flow channels during the exploitation of natural gas hydrates.

In this study, a comprehensive gas hydrates exploitation model was established. This model takes into account various critical factors, including ice formation, the migration of hydrate and ice particles, and the potential occurrence of flow channel jamming. The influence of solid particles migration and jamming on the gas production and the underlying mechanisms was revealed.

Models

In this study, the Two-fluid Model (TFM) and Finite Volume Method (FVM) were used to simulate the process of hydrate exploitation. The gas hydrates exploitation model can be established with the ice formation, solid migration, and channel jamming considered. Assuming that the skeleton of porous media does not deform and the sand is stationary, there are four phases (water, gas, hydrate, and ice) migrations in this model. The convection terms associated with hydrate and ice phases are added to the conservation equations of mass and energy.

Mass conservation equation:

$$-\mathrm{\nabla }\cdot (\varphi S_w{\rho }_w{\overrightarrow{v}}_w)-\mathrm{\nabla }\cdot (\varphi S_g{\rho }_g{\overrightarrow{v}}_g)-\mathrm{\nabla }\cdot (\varphi S_h{\rho }_h{\overrightarrow{v}}_{hi})-\mathrm{\nabla }\cdot (\varphi S_i{\rho }_i{\overrightarrow{v}}_{hi})+{\dot{m}}_w+{\dot{m}}_g-{\dot{m}}_h={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}(\varphi S_w{\rho }_w+\varphi S_g{\rho }_g+\varphi S_h{\rho }_h+\varphi S_i{\rho }_i)}$$

(1)

Energy conservation equation:

$$\begin{array}{rl}& \mathrm{\nabla }\cdot \{[(1-\varphi ){\lambda }_s+\varphi S_w{\lambda }_w+\varphi S_g{\lambda }_g+\varphi S_h{\lambda }_h+\varphi S_i{\lambda }_i]\mathrm{\nabla }T\}-{\dot{m}}_h\mathrm{\Delta }{H}_h+{\dot{m}}_i\mathrm{\Delta }{H}_i\\ & -\mathrm{\nabla }\cdot (\varphi S_w{\rho }_w{c}_w{\overrightarrow{v}}_{w}T)-\mathrm{\nabla }\cdot (\varphi S_g{\rho }_g{c}_g{\overrightarrow{v}}_{g}T)-\mathrm{\nabla }\cdot (\varphi S_h{\rho }_h{c}_h{\overrightarrow{v}}_{hi}T)-\mathrm{\nabla }\cdot (\varphi S_i{\rho }_i{c}_i{\overrightarrow{v}}_{hi}T)\\ & ={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}\{[(1-\varphi ){\rho }_s{c}_s+\varphi S_w{\rho }_w{c}_w+\varphi S_g{\rho }_g{c}_g+\varphi S_h{\rho }_h{c}_h+\varphi S_i{\rho }_i{c}_i]T\}}\end{array}$$

(2)

The fluid velocity is calculated by Darcy's law. However, the solid velocity is obtained by solving the momentum equation:

$$-\mathrm{\nabla }\cdot [\varphi (S_h+S_i){\rho }_{hi}{\overrightarrow{v}}_{hi}{\overrightarrow{v}}_{hi}]-\varphi (S_h+S_i)\mathrm{\nabla }{P}_{hi}+\mathrm{\nabla }\cdot \stackrel{-}{\stackrel{-}{\tau }}+\varphi (S_h+S_i){\rho }_{hi}\overrightarrow{g}+{\beta }_{ls}({\overrightarrow{v}}_{wg}-{\overrightarrow{v}}_{hi})-{\dot{m}}_h{\overrightarrow{v}}_{hi}+{\dot{m}}_i{\overrightarrow{v}}_{hi}={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}[\varphi (S_h+S_i){\rho }_{hi}{\overrightarrow{v}}_{hi}]}$$

(3)

where $\stackrel{=}{\tau }$ is the solid phase shear stress, which can be calculated by:

$$\stackrel{=}{\tau }=[-P_{hi}+\varphi (S_h+S_i){\lambda }_{hi}\mathrm{\nabla }\cdot {\overrightarrow{v}}_{hi}]\stackrel{-}{I}+\varphi (S_h+S_i){\mu }_{hi}\left[\mathrm{\nabla }{\overrightarrow{v}}_{hi}+{(\mathrm{\nabla }{\overrightarrow{v}}_{hi})}^T-{\displaystyle \frac{2}{3}(\mathrm{\nabla }\cdot {\overrightarrow{v}}_{hi})\stackrel{-}{I}}\right]$$

(4)

where ${\lambda }_{hi}$ is the solid particle bulk viscosity, which can be calculated by:

$${\lambda }_{hi}={\displaystyle \frac{4}{3}\varphi (S_h+S_i){\rho }_{hi}{d}_{hi}{g}_{hi}(1+e_{hi}){\left({\displaystyle \frac{{\mathrm{\Theta }}_{hi}}{\pi }}\right)}^{\frac{1}{2}}}$$

(5)

where $e_{hi}$ is the collision coefficient of restitution of solid particles (assuming that the solid particles are elastic collision, $e_{hi}=1.0$ in this study); $g_{hi}$ is the radial distribution function of solid particles:

$$g_{hi}={\left[1-\varphi \left({\displaystyle \frac{S_h+S_i}{{\alpha }_{max}}}\right)\right]}^{-1}$$

(6)

where the value of maximum packing limit of solid particles ${\alpha }_{max}$ is 0.74 (Chen, 2014).

In equation(4), the shear viscosity of the solid phase, ${\mu }_{hi}$ , consists of the collision term and the dynamic term (the friction term between particles is ignored).

$${\mu }_{hi}={\mu }_{hi,col}+{\mu }_{hi,kin}$$

(7)

$${\mu }_{hi,col}={\displaystyle \frac{3}{5}{\lambda }_{hi}}$$

(8)

$${\mu }_{hi,kin}={\displaystyle \frac{\varphi (S_h+S_i){\rho }_{hi}{d}_{hi}\sqrt{{\mathrm{\Theta }}_{hi}\pi }}{6(3-e_{hi})}\left[1+{\displaystyle \frac{2}{5}\varphi (S_h+S_i)(1+e_{hi})(3e_{hi}-1)g_{hi}}\right]}$$

(9)

For the solid phase pressure, $P_{hi}$ , it is measured by particle pseudothermal temperature, ${\mathrm{\Theta }}_{hi}$ :

$$P_{hi}=\varphi (S_h+S_i){\rho }_{hi}{\mathrm{\Theta }}_{hi}+2{\rho }_{hi}(1+e_{hi}){\varphi }^2(S_h+S_i{)}^2{g}_{hi}{\mathrm{\Theta }}_{hi}$$

(10)

$$\begin{array}{rl}& {\displaystyle \frac{3}{2}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[\varphi \left(S_h+S_i\right){\rho }_{hi}{\mathrm{\Theta }}_{hi}\right]+{\displaystyle \frac{3}{2}\mathrm{\nabla }\cdot \left[\varphi \left(S_h+S_i\right){\rho }_{hi}{\overrightarrow{v}}_{hi}{\mathrm{\Theta }}_{hi}\right]}}}\\ & =\overline{\overline{\tau }}:\mathrm{\nabla }{\overrightarrow{v}}_{hi}+\mathrm{\nabla }\cdot \left(k_{\mathrm{\Theta }}\mathrm{\nabla }{\mathrm{\Theta }}_{hi}\right)-3{\beta }_{ls}{\mathrm{\Theta }}_{hi}\end{array}$$

(11)

where $k_{\mathrm{\Theta }}$ is the transport coefficient of the solid phase pseudothermal temperature, which can be calculated by:

$$k_{\mathrm{\Theta }}={\displaystyle \frac{150{\rho }_{hi}{d}_{hi}\sqrt{{\mathrm{\Theta }}_{hi}\pi }}{384(1+e_{hi})g_{hi}}{\left[1+{\displaystyle \frac{6}{5}\varphi (S_h+S_i)g_{hi}(1+e_{hi})}\right]}^2+2{\rho }_{hi}{d}_{hi}{g}_{hi}{\varphi }^2(S_h+S_i{)}^2\sqrt{{\displaystyle \frac{{\mathrm{\Theta }}_{hi}}{\pi }}}}$$

(12)

The momentum exchange coefficient between solid and fluid phases, ${\beta }_{ls}$ , is calculated by the Syamlal-O’Brien model:

$${\beta }_{ls}={\displaystyle \frac{3{\varphi }^2(S_h+S_i)(S_w+S_g){\rho }_{wg}}{4v_{r,hi}^2{d}_{hi}}{C}_D\left({\displaystyle \frac{R{e}_{hi}}{v_{r,hi}}}\right)|{\overrightarrow{v}}_{hi}-{\overrightarrow{v}}_{wg}|}$$

(13)

$$v_{r,hi}=0.5\left[A-0.06{\mathrm{Re}}_{hi}+\sqrt{{(0.06R{e}_{hi})}^2+0.12{\mathrm{Re}}_{hi}(2B-A)+A^2}\right]$$

(14)

where the solid Reynolds number, $R{e}_{hi}$ , and the coefficient of drag force, $C_D$ , are calculated as the following formulas:

$${\mathrm{Re}}_{hi}={\displaystyle \frac{{\rho }_{wg}{d}_{hi}(S_w+S_g)|{\overrightarrow{v}}_{hi}-{\overrightarrow{v}}_{wg}|}{S_w{\mu }_w+S_g{\mu }_g}}$$

(15)

$$C_D={\left(0.63+4.8\sqrt{{\displaystyle \frac{v_{r,hi}}{R{e}_{hi}}}}\right)}^2$$

(16)

In order to solve the solid phase momentum equation to get the solid phase velocity, it needs to assume a velocity value of the solid phase to calculate the coefficients in equation (3). Then, the modified value of the solid phase velocity is obtained. And then, through iterative calculation, the velocity of the solid phase in the current time step is determined when the convergence condition is satisfied. At the same time, there is another criterion for the solid velocity. When jamming occurs within the flow channels of the computational grid, the solid phase velocity at that location is zero. Moreover, this criterion has the highest priority in controlling the velocity values. Whether jamming occurs in the flow channels is determined by the transient jamming rate calculation formula derived by Lafond et al. (Lafond, 2014), as shown in equation (17). When the jamming rate goes over 1, the flow channel in the control volume is determined to be jamming, and the solid phase velocity is forced to be 0 in the model.

$$r(t)=(-\dot{N}\ln p{)}_t=-{\displaystyle \frac{3{\phi }_C{v}_f{d}_{pipe}^2{(R_r-1)}^2\ln (1-e^{-\alpha (R_r^2-1)})}{2d_o^2{d}_p}}$$

(17)

During the exploitation of NGH, a solid particle bed is formed within the flow channels due to the effects of solid–liquid density differences and adhesion forces (Duan et al., 2024). Some idealized assumptions are included to apply equation (17) to the numerical simulation of hydrate exploitation: (1) the flow channels in porous media are simplified as two-dimensional rectangles; (2) there is one and only one flow channel in each control volume; (3) upon hydrate dissociation, ice formation and melting, only the number of particles change, while the size of solid particles remained the same, and the particles are equal diameter spheres; (4) the diameter of solid particles is taken as 50 µm (Grasso, 2015); (5) in each control volume with solid phase, solid particles are completely deposit, that is, the length of the particle bed L is equal to the length of the grid $L_{net}$ , that is, $L=L_{net}$ ; (6) the hydrate in porous media is filling type. Based on these assumptions, Figure 1 shows the schematic diagram of flow channel and particle bed dimensions in a two-dimensional single control volume. According to the assumptions (1) and (2), the diameter of flow channel ${\mathrm{d}}_{pipe}$ remains consistent in the uniform grid, ${\mathrm{d}}_{pipe}=\varphi L_{net}$ . The orifice size $d_o$ formed by the particle bed is given by:

$$d_o=\varphi L_{net}-{\displaystyle \frac{2\varphi (S_h+S_i)L_{net}^2}{L}}$$

(18)

Figure 1.

Schematic diagram of flow channel and particle bed size in one control volume under idealized assumptions.

The coupled iterative process is illustrated in the computational flowchart as Figure 2. Therefore, the interaction between solid phase velocity, solid phase saturation, and jamming rate is constructed to describe the influence of solid phase displacement on hydrate exploitation process. Using this model, the influence of solid phase migration and channel jamming on the gas production from hydrate exploitation is studied.

Figure 2.

Computational flowchart of the coupled iterative process.

At the beginning of every time step, the parameters obtained in the previous step are substituted into equation (17) to compute $r(t)$ . If $r(t)\ge 1$ , the solid phase velocity ${\overrightarrow{v}}_{hi}$ is set to zero. Equation (1) is then discretized on a uniform grid. Darcy's law is inserted to replace fluid velocity by pressure and permeability, yielding a system of linear equations for pressure. Phase properties, saturations, and permeabilities evaluated at the previous time step are used to assemble the coefficients, and the pressure field is obtained iteratively. Phase saturations are subsequently updated through the hydrate-decomposition kinetics model, the Stefan scheme, and the individual mass-balance equations. Finally, the temperature field is solved from the linear system produced by discretizing equation (2). If $r(t)<1$ , the solid velocity from the previous time step is first inserted into equation (3), and a corrected solid velocity is obtained by discrete iteration. This corrected velocity is then fed back into equation (1) and equation (2) to re-solve the pressure and temperature fields. Upwind differencing is employed for the pressure and temperature equations. The tridiagonal matrix algorithm is adopted to solve the discretized linear equation system in a fully implicit format. A fixed pressure is imposed at the production well, while fixed-temperature conditions are prescribed on the remaining three boundaries.

Results and analysis

Model validation

This study validated the numerical model using experimental data obtained from marine sediment drilled from the South China Sea. The initial conditions of the numerical simulation were configured to match the experimental conditions reported in Ref. (Wang et al., 2023), as detailed in Table 1. Also, we compared the numerical simulation results with the experimental data of the prepared core samples, with the initial parameter settings shown in Table 2.

Table 1.

Initial conditions of the Wang’s experiments (Wang et al., 2023).

Parameter	Value	Parameter	Value
Reservoir size (mm × mm)	38 × 76	Initial hydrate saturation, $S_h$ (%)	22.96
Initial permeability, $K_0$ (md)	18	Initial water saturation, $S_w$ (%)	30
Porosity, ϕ (%)	42	Initial gas saturation, $S_g$ (%)	47.04
Initial reservoir pressure, $P_0$ (MPa)	8	Outlet pressure, $P_W$ (MPa)	3
Initial reservoir temperature, $T_0$ (K)	277.15	Permeability reduction index, $N$	4

Table 2.

Initial conditions of the Tang’s experiments (Tang et al., 2007).

Parameter	Value	Parameter	Value
Reservoir size (cm × cm)	38 × 38	Initial hydrate saturation, $S_h$ (%)	17.6
Initial permeability, $K_0$ (md)	1970	Initial water saturation, $S_w$ (%)	38.9
Porosity, ϕ (%)	40	Initial gas saturation, $S_g$ (%)	43.5
Initial reservoir pressure, $P_0$ (MPa)	3.24	Outlet pressure, $P_W$ (MPa)	2.25
Initial reservoir temperature, $T_0$ (K)	274.85	Permeability reduction index, $N$	4

The simulation results of the cumulative gas production with and without solid phase migration are compared with the experimental results. As shown in Figure 3, when considering solid phase migration in the model, the gas production rate significantly decreased. When the solid phase is assumed to be stationary in the numerical simulation, the average relative error of the cumulative gas production compared with the experimental data is 14.3%. When the TFM that accounts for the migration of hydrate and ice is employed, this error is reduced to 6.2%. The discrepancy in the numerical result is mainly attributed to the fact that the TFM model accounts only for the migration of hydrate and ice, whereas the movement of sand particles and the deformation of the porous-medium skeleton are neglected (Figure 4).

Figure 3.

Comparison of the cumulative gas production using the TFM, the model with solid stationary and experimental data (Wang et al., 2023).

Figure 4.

Comparison of the reservoir gas rate using the TFM, the model with solid stationary and experimental data (Tang et al., 2007).

The simulation results of the gas production rate with and without solid phase migration are compared with the experimental results. As shown in Figure 3, when considering solid phase migration in the model, the gas production rate significantly decreased. Upon analyzing the degree of fit between the calculated results of the two models and the experimental data, the correlation coefficient R2 of the green curve (representing the gas production rate curve considering solid phase migration) is 0.939, and that of the blue curve (representing the gas production rate curve considering solid phase stationary) is 0.681. Above all, it suggests that the model considering solid phase migration describes the hydrate exploitation and predicts the gas production better.

Influence of solid phase migration on reservoir absolute permeability

Previous studies have indicated that the key factors influencing the gas production rate are the dissociation rate of hydrate and the capacity of gas and liquid in the reservoir to flow through the porous media into the production well (Konno et al., 2010). These two factors are mainly governed by permeability, temperature, and pressure in the reservoir. Below, numerical calculations are performed using the parameter values in Table 2 as initial conditions to discuss the influence of solid phase migration on these aspects.

The absolute permeability of the reservoir is a function of the saturation of hydrate and ice phases. Changes in saturation are affected by phase changes and solid migration. Experimental data indicate that the reservoir's permeability decreases with increasing fluid velocity (Gruesbeck and Collins, 1982). The reason is that the larger fluid velocity increases the drag force on the solid particles, thus strengthening the migration effect of the solid particles (Syamlal and O’Brien, 1989). Comparing the permeability changes with and without solid phase migration, as shown in Figure 5, the permeability increases faster in the solid stationary model, especially in the area close to the mining well (the left edge of the computational domain). According to Darcy's law for flow, fluid velocity rises with an increase in pressure gradient. The pressure of the producing well leads to a larger pressure gradient near the production well. Considering the migration of hydrate phase and ice phase, the solid phase saturation is high at the same time, and the reservoir permeability is low. Consequently, there is a significant difference in permeability between the two models in the vicinity of the well. This is consistent with the phenomenon observed by Gruesbeck et al. (Gruesbeck and Collins, 1982) in their experiments. Based on the comparison results, the lower permeability in the near-wellbore area caused by the migration of hydrate and ice phases is one of the reasons for gas production inhibition.

Figure 5.

Changes of the calculated reservoir absolute permeability over time using the model with solid stationary and solid migration (results from the model with solid stationary: (a) 5 min, (b) 10 min, (c) 15 min; results from the model with solid migration: (d) 5 min, (e) 10 min, (f) 15 min).

At the same time, the poor gas–liquid flow caused by low permeability will also be applied to the temperature field variation of the reservoir through convective heat transfer, thus affecting the dissociation rate of hydrate. In turn, it inhibits the permeability growth.

Influence of solid phase migration on temperature and pressure

The influence of solid phase migration on temperature and pressure was also studied. To monitor temperature and pressure changes at different locations in the reservoir, four monitoring points were placed in the calculation domain as shown in Figure 6. The pressure at point 1, located on the wall of the mining well, remains consistent with the well pressure. Therefore, it makes sense to monitor the pressure over time at points 2 and 4.

Figure 6.

Diagram of the location of the monitoring points in the computational domain.

Figure 7 shows that the pressure curves at the same points exhibit close alignment between the models with and without solid phase migration, especially in the early stages of production. As the production proceeds, the pressure obtained by the model with solid phase migration is slightly higher. This is due to the small increase in local pressure caused by the decrease of gas and liquid flow in the reservoir due to the accumulation of the solid phase near the production well area. The increase in local pressure tends to inhibit the decomposition rate of hydrate.

Figure 7.

Comparison of the pressure change at points 2 and 4 using the model with and without solid migration.

The changes in the reservoir's temperature field calculated by the models with and without solid migration are shown in Figure 8. Based on Darcy's law for flow through porous media, under the same pressure gradient, the migration of solid particles results in a lower absolute permeability of the reservoir and a lower velocity of gas–liquid flow. In the energy conservation equation (7), the decreased fluid velocity leads to a reduction in convective heat transfer between the control volumes, ultimately causing alterations in the calculated reservoir temperature field. The changes in the temperature field are primarily controlled by convective heat transfer induced by gas–liquid flow and the endothermic dissociation of hydrates. Therefore, this model uses an upwind scheme for numerical computations. When the flow channel is unobstructed (without considering the migration of solid-phase particles), enthalpy flow from upstream is transmitted downstream. As the region near the extraction well experiences lower pressure and more hydrate dissociation heat, the two effects counterbalance each other, resulting in a relatively uniform temperature field. However, when considering the migration of solid-phase particles, these particles accumulate near the extraction well, reducing the gas–liquid flow capacity in that region and decreasing the convective heat transfer rate. The heat absorbed by hydrate dissociation in that area cannot be promptly replenished by the incoming energy from upstream, leading to a faster temperature drop. Simultaneously, enthalpy flow converges in the central region of the reservoir, causing higher temperatures in that area. The early formation of low temperature inhibits the dissociation rate of hydrate near the extraction well, which, in turn, leads to a reduction in the gas production rate of the reservoir.

Figure 8.

Changes of the calculated temperature over time using the model with and without solid migration (results from the model without solid migration: (a) 5 min, (b) 10 min, (c) 15 min; results from the model with solid migration: (d) 5 min, (e) 10 min, (f) 15 min).

At the same time, since the initial reservoir temperature was set near the freezing point, the drop in temperature resulted in the formation of the ice phase. The formation of the ice phase is mainly controlled by temperature variation. Consequently, in the model considering solid migration, the ice forms relatively earlier in the area near the mining well as shown in Figure 9. Unlike the calculation results of the model with solid stationary, when the solid migration is taken into account, the ice phase becomes more concentrated in the near-wellbore area, with higher ice phase saturation is higher concurrently. On one hand, the temperature in this region drops rapidly and the amount of ice is more. On the other hand, the ice phase generated slightly further away from the mining well will move toward it, resulting in increased ice phase saturation in this area. This is one of the reasons for the decrease in permeability in the area close to the well.

Figure 9.

Changes of the calculated ice distribution over time using the model with and without solid migration (results from the model without solid migration: (a) 5 min, (b) 10 min, (c) 15 min; results from the model with solid migration: (d) 5 min, (e) 10 min, (f) 15 min).

Influence of jamming on reservoir absolute permeability

For the process of gas hydrate production by depressurization, lower well pressure increases the rate of hydrate dissociation in the reservoir, particularly in the vicinity of the mining well (Yu et al., 2018a). This should lead to a reduction in solid phase saturation and an elevation in reservoir absolute permeability. However, when accounting for solid migration and flow jamming, the permeability pattern is different.

In this study, the reservoir absolute permeability at monitoring points 1 and 3 was calculated under different mining well pressure (2.25 MPa and 1.5 MPa). As shown in Figure 10, at the same reservoir location, the permeability is higher with lower well pressure. This is due to the factor that lower well pressure promotes hydrate dissociation, resulting in a greater decrease in total solid saturation. However, under a well pressure of 1.5 MPa, the permeability curve at monitoring point 1 has a special variation trend. Initially, the permeability experiences a rapid increase, followed by gradual steady growth, and then a sudden decrease, dropping even below the initial permeability value. The rapid initial increase in permeability is caused by the extensive dissociation of hydrates at the boundary. Subsequently, the solid phase continues to migrate toward the production well area under a large pressure gradient, leading to a gradual slowdown in the rate of permeability growth at point 1. Along with the continuous accumulation of solid phase, the jamming judgment condition is triggered. The jamming of the flow channel halts the solid phase there, while the movement of the solid phase from adjacent areas rapidly increases the solid phase saturation, resulting in a sudden and significant decrease in permeability. With the further dissociation of hydrate, the jamming of the flow channel is soon removed and the permeability rises again. The occurrence of this phenomenon indicates that the low well pressure condition is easy to promote near-well area jamming. The jamming leads to an abnormal reduction in permeability in the area, contributing to the decline in the reservoir's gas production rate. The pressure field inside the reservoir also changes accordingly. As shown in Figure 11, when the mining well pressure is set at 1.5 MPa, the rate of pressure decline at monitoring point 2 (located near the mining well) slows around 30 min and then accelerates at approximately 50 min. Consequently, within this time interval, the pressure difference between the regions at point 3 and point 2 is further reduced, diminishing the gas–liquid flow driving force. Therefore, it is suggested that appropriate well pressure should be set based on reservoir conditions during hydrate depressurization production. Excessive reduction in well pressure could result in decreased production efficiency.

Figure 10.

Comparison of the reservoir absolute permeability at points 1 and 3 under different mining well pressure (2.25 MPa and 1.5 MPa).

Figure 11.

Comparison of the pressure change at points 1 to 4 using the model with solid migration under the mining well pressure 1.5 MPa.

Inhibition mechanism on the gas production rate

According to the aforementioned discussion, the migration of solid particles and flow channel jamming significantly influence the variations in key parameters during hydrate exploitation. Changes in major parameters such as permeability, temperature, and pressure ultimately manifest in the decomposition gas production rate of hydrates. As depicted in Figure 12, the gas production rate calculated by this model initially experiences a rapid decrease, followed by gradual increases with minor fluctuations, eventually stabilizing. Notably, throughout the process, the gas production rates obtained from this model are consistently lower than those from the model that does not account for solid particle migration and flow channel jamming.

Figure 12.

Comparison of the gas production rate between models considering (TFM) and not considering solid migration and flow channel jamming.

Evidently, the migration of solid particles and flow channel jamming exert a notable inhibitory effect on hydrate decomposition and gas production. This inhibitory mechanism aligns with the earlier analysis of variations in key parameter characteristics. First, solid particle migration leads to lower permeability in the near-wellbore region, reducing fluid output rates and suppressing reservoir gas production rates. Moreover, excessively low wellbore pressures can intensify solid particle migration, potentially inducing flow channel jamming, which results in periodic and substantial drops in reservoir permeability, thereby enhancing the volatility of gas production rates. Second, the accumulation of solid particles due to migration in the near-wellbore area diminishes the flow capacity for gas–liquid movement, decreasing convective heat exchange. As a consequence, the heat absorbed during hydrate decomposition is not adequately replenished, leading to temperature reductions in that region. Subsequently, this decrease in temperature leads to a decline in the rate of hydrate decomposition and, consequently, in gas production rates. Simultaneously, when local temperatures drop below the freezing point, ice formation occurs, further diminishing the rates of hydrate decomposition and reservoir gas production. Third, the decrease in gas–liquid flow capacity caused by solid particle migration results in a minor increase in internal reservoir pressure, somewhat inhibiting the rate of hydrate decomposition. However, the reduction in hydrate decomposition gas production counteracts the rise in pressure, leading to relatively small changes in reservoir pressure and limited inhibitory effects on gas production rates.

Conclusion

In this study, a comprehensive model of hydrate exploitation was established, taking into account solid migration and channel jamming. This model reasonably describes heat and mass transfer among the four components (hydrate, gas, water, and ice) and the three phases (solid, liquid, and gas). By considering solid phase migration, the calculated results of the model show a decrease in the gas production rate, and the consistency of the gas production rate between numerical simulation and experiments was significantly improved. Through this model, the mechanism of the inhibition effect of solid migration and ice formation on the gas production was revealed. First, the influence of hydrate phase and ice phase migration and jamming inhibits the absolute permeability of the reservoir due to the accumulation of solid particles near the well area. Second, due to the jamming caused by solid particles, fluid flow is inhibited, and convective heat transfer is reduced, resulting in a rapid decline in local reservoir temperature. Consequently, near the well area, the temperature is more likely to drop below freezing, leading to an ice phase formation. Due to solid phase migration and flow channel jamming, the aforementioned parameters undergo changes, ultimately resulting in a decrease in the gas production rate and gas hydrate dissociation rate. Therefore, measures should be taken to suppress the migration of the hydrate phase in the process of hydrate exploitation to enhance the gas production rate of the reservoir.